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Theoretical  Background 

1.1  Introduction 

A  method  for  computing  indices  of  muscle  fatigue  is  introduced  by  decomposing  the  surface 
electromyographic  (SEMG)  signal  using  a  filter  bank  and  estimating  from  the  subbands  of  the 
filter  bank  both  by  the  instantaneous  frequency  (IF)  and  the  instantaneous  amplitude  (IA).  These 
time  dependent  parameters  are  then  used  to  compute  indices  of  muscle  increase  force,  recovery, 
decrease  force  and  fatigue. 

Surface  electromyographic  signals  have  been  proposed  as  a  potential  tool  to  detect  and  measure 
muscle  fatigue  [8],  [11],  [13],  [3].  Specifically,  muscle  fatigue  involves  a  decrease  in  the 
maximum  force  and  a  measurable  shift  in  the  power  spectrum  of  the  SEMG  signal,  which  is 
commonly  reported  to  occur  toward  lower  frequencies  [15],  [20],  [16],  although  there  are  fewer 
reports  of  increased  frequencies  [16].  Consequently,  the  mean  frequency,  median  frequency,  and 
IF  have  been  extensively  used  in  the  assessment  of  muscle  fatigue  as  indicators  of  how  the 
frequency  of  the  SEMG  signal  changes  over  time  [20],  [25],  [30],  [19].  Additionally,  muscle 
fatigue  also  involves  an  increase  in  the  SEMG  signal’s  amplitude  [15],  [20],  [32],  [9].  The  root 
mean  square,  the  average  rectified  value,  and  a  recent  technique  introduced  by  Merletti  [21],  [15], 
[19],  [27]  have  been  proposed  as  amplitude  estimators  of  the  SEMG  signal,  that  can  then  be  used 
as  an  additional  input  to  the  analysis  of  muscle  fatigue. 

When  considering  changes  in  the  frequency  and  amplitude  of  the  SEMG  signal  to  investigate 
muscle  fatigue,  the  level  of  stationarity  of  the  signal  should  be  considered.  Different  analysis 
techniques  have  been  found  appropriate  for  assessment  of  muscle  fatigue,  but  they  are  dependent 
on  the  level  of  stationarity  [9],  [12],  [25],  [11].  The  SEMG  signal  is  considered  wide-sense 
stationary  when  recorded  during  isometric  constant  force  contractions  over  time  intervals  lasting 
from  one  to  two  seconds  [21].  In  contrast  the  SEMG  signal  is  considered  non-stationary  when 
recorded  during  dynamic  contractions.  Non-stationarities  of  the  SEMG  signal  can  be  classified  as 
slow  or  fast  [4] . 

Time-frequency  (TF)  distributions  from  the  Cohen  class  and  wavelets  have  demonstrated  to  be 
useful  approaches  in  the  analysis  of  SEMG  signals  with  static  and  dynamic  muscle  contractions 
[8],  [12],  [25],  [11],  [4],  [13],  [3].  The  spectrogram  is  an  example  of  TF  distributions  that  has 
been  used  in  the  analysis  of  static  and  dynamic  muscle  contractions  with  low  to  medium 
non-stationarities  [32],  [9].  The  objective  of  TF  distributions  and  wavelets  in  this  context  is  to 
provide  information  about  how  the  energy  of  the  signal  is  distributed  over  time  and  frequency.  In 
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most  cases,  they  have  been  used  to  estimate  the  frequency  but  not  the  amplitude.  An  exception  is 
the  work  by  Luttman  [20],  [19],  which  introduces  a  joint  analysis  of  spectrum  and  amplitude 
(JASA)  of  SEMG  signals.  During  this  research,  JASA  is  referred  as  a  joint  analysis  of  frequency 
and  amplitude. 

In  order  to  assess  muscle  fatigue,  the  IF  and  the  IA  from  the  subband  components  of  SEMG 
signals  are  estimated  and  decomposed  with  a  32-channel  cosine  modulated  filter  bank.  The 
SEMG  signals  provided  by  the  Air  Force  Research  Laboratory  (AFRL)  were  recorded  from  26 
healthy  human  subjects  during  prolonged  wear  of  weighted  flight  helmets.  A  linear  regression 
model  is  adopted  to  compute  the  slopes  from  the  IF  and  IA  estimates.  The  slopes  derived  from  the 
filter  bank  are  then  compared  with  the  slopes  derived  from  the  spectrogram  and  the  smoothed 
pseudo  Wigner-Ville  distribution  (SPWVD).  Finally,  the  slopes  are  used  as  indices  of  muscle 
increase  force,  recovery,  decrease  force  and  fatigue,  and  are  correlated  to  perceived  levels  of 
discomfort  reported  by  the  subjects. 

1.2  Overview 

An  extensive  search  of  the  literature  shows  that,  according  to  most  researchers,  in  muscle  fatigue 
there  is  a  shift  of  the  total  power  spectrum  of  the  SEMG  signal  toward  lower  frequencies  because 
the  conduction  velocity  of  the  motor  actions  potentials  on  the  muscle  membrane  decreases. 
Spectral  variables  such  as  the  mean  frequency  (MNF),  the  median  frequency  (MF)  and  the 
instantaneous  frequency  (IF)  have  been  extensively  proposed  to  measure  this  shift  in  the  power 
spectrum.  However,  it  is  still  unclear  which  spectral  variable  is  the  best  in  terms  of  detecting  this 
shift.  For  instance  some  authors  state  that  the  MF  is  preferred  since  it  is  less  sensitive  to  noise. 
Rather,  others  state  that  the  MNF  is  more  stable  and  sensitive  to  changes  in  the  underlying 
spectrum  [20].  Also,  it  has  been  proposed  that  in  the  assessment  of  muscle  fatigue  one  should 
analyze  changes  in  both  the  power  spectrum  and  the  amplitude  of  the  SEMG  signal  [15],  [20], 
[16],  [19],  [32]  [9], 

In  the  literature  different  techniques  are  found  to  estimate  spectral  variables  (MNF,  MF  and  IF) 
and  the  amplitude  of  the  SEMG  signal.  A  classical  method  for  estimating  the  MNF  and  the  MF  is 
explained  in  [15],  [20],  [19],  [32],  [9].  First,  the  SEMG  signal  is  divided  into  equal  segments 
(epochs)  with  a  duration  of  500  ms  or  1  s.  If  the  signal  is  non-stationary,  then  the  SEMG  signal  is 
divided  into  epochs  of  250  ms  or  500  ms  [9].  On  the  other  hand,  if  the  signal  is  highly 
non-stationary,  then  the  signal  is  divided  into  epochs  shorter  than  125  ms.  Some  authors  suggest 
an  overlap  in  the  epochs  of  50%  to  reduce  the  variance  [20] .  However,  epoch  overlapping  is  not 
recommended  since  it  increases  the  computational  time  without  providing  significant 
improvement  of  the  quality  of  the  estimates  [9]. 

After  the  SEMG  signal  is  divided  into  different  epochs,  the  power  spectral  density  (PSD)  of  the 
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epoch  is  then  estimated  by  applying  the  absolute  value  squared  of  the  fast  Fourier  transform 
(FFT).  The  MNF  and  the  MF  of  a  SEMG  signal  can  be  estimated  as 


MNF 


£°oof(S(f))2df 
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1  f'00 

MF  =  ~  ( S(f))2df ,  (1.2) 
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where  S(f)  is  the  estimated  PSD  of  the  SEMG  signal  s(t),  and  /  is  the  frequency  variable. 
Finally,  the  amplitude  of  the  SEMG  signal  can  be  estimated  from  each  epoch  as 
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where  T  is  the  length  of  the  interval  where  the  SEMG  signal  s(t)  is  analyzed  [15],  [19],  [32]  [9]. 
Another  technique  in  the  literature  regarding  the  estimation  of  spectral  variables  (MNF,  MF,  and 
amplitude)  is  the  use  of  time- varying  autoregressive  models  (TVAR).  SEMG  signals  can  be 
represented  using  an  autoregressive  (AR)  model  with  parameters  changing  with  time.  Thus,  the 
time-varying  spectrum  can  be  estimated  from  the  time-varying  model  parameters,  and  the  IF  of 
the  non- stationary  signal  can  be  extracted  [31],  [33].  Furthermore,  the  selection  of  the  model 
order  for  an  AR  is  not  critical,  and  an  order  of  10  seems  to  be  appropriate  for  epochs  of  125ms, 
250ms,  500ms,  and  Is  [9].  TVAR  can  be  used  to  monitor  the  muscle  fatigue  development  during 
isometric  contraction  [33],  [31]. 

Some  authors  propose  to  estimate  the  IF  by  computing  the  first  derivative  of  the  phase  of  the 
Hilbert  transform  of  a  signal  [10].  However,  other  authors  point  out  that  the  IF  derived  from  this 
method  is  meaningful  just  for  monocomponent  signals.  In  the  case  of  multicomponent  signals 
such  as  the  SEMG  signals,  one  would  obtain  at  some  points  negative  frequencies  and  frequencies 
outside  the  bandwidth  of  the  signal  [30].  As  a  result  this  method  is  not  suggested  to  estimate  the 
IF  of  SEMG  signals.  An  alternative  solution  to  this  issue  is  proposed  in  [30]  and  [10].  In  [30],  a 
technique  known  as  the  Hilbert-Huang  transform  (HHT)  is  proposed.  The  HHT  is  based  on 
decomposing  the  SEMG  signal  into  a  set  of  intrinsic  mode  functions  in  such  a  way  that  they  may 
be  Hilbert  transformed.  It  is  concluded  that  the  HHT  is  the  best  choice  in  the  analysis  of  muscle 
fatigue,  and  that  wavelets  are  lowest  in  performance  after  a  comparison  of  the  slopes  derived  from 
HHT,  TVAR,  and  wavelets.  The  other  solution  given  in  [10]  proposes  to  extract  the  IF  from  an 
adaptive  TF  distribution  of  the  signal.  This  adaptive  TF  distribution  is  obtained  by  decomposing 
the  signal  into  components  with  good  TF  localizations  and  combining  the  Wigner  distribution  of 
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the  components.  Hence,  the  proposed  TF  distribution  is  free  of  cross  terms,  positive,  and  with 
good  TF  localization. 

Recently,  TF  analysis  has  been  adopted  in  the  analysis  of  both  stationary  and  non-stationary 
electromyographic  signals.  The  fact  that  spectral  variables  can  be  extracted  from  a  TF  distribution 
makes  this  an  interesting  distribution  in  the  assessment  of  muscle  fatigue.  In  the  literature 
different  approaches  to  estimate  a  TF  distribution  is  found.  Included  among  them  are  the 
Wigner-Ville  distribution  (WVD),  the  short-time  Fourier  transform  (STFT)  or  spectrogram,  the 
Choi- Williams  distribution  (CWD),  the  smoothed  pseudo  Wigner-Ville  distribution  (SPWVD), 
and  wavelets.  A  comparison  of  the  WVD,  spectrogram,  SPWVD  and  CWD  is  attained  in  [4],  [3], 
[8]  and  [1].  The  results  show  that  the  WVD  is  an  efficient  tool  in  the  analysis  of  monocomponent 
signals;  however,  its  performance  in  the  analysis  of  multi-component  signals  is  negatively 
affected  by  cross  terms  in  the  estimates  of  the  spectral  variables.  The  spectrogram  is  restricted  by 
the  fact  that  a  good  localization  in  time  leads  to  a  poor  localization  in  frequency  and  vice-versa; 
furthermore,  it  assumes  the  analyzed  signal  to  be  stationary  inside  each  considered  fixed  time 
window,  which  is  not  true  in  the  case  of  dynamic  contractions.  Regarding  the  SPWVD,  it  destroys 
components  of  the  signal  in  the  rejection  of  the  cross  terms  and  also  employs  fixed  time  and 
frequency  windows.  Finally,  the  CWD  seems  to  be  very  effective  in  decreasing  the  effects  of 
cross  terms  and  in  retaining  most  of  the  useful  properties  of  a  TF  distribution;  however,  it  does  not 
dramatically  improve  the  results  with  respect  to  the  SPWVD  in  terms  of  better  peak  tracking 
capabilities. 

Wavelet  and  Wavelet  packets  (WP)  have  been  recently  introduced  in  the  analysis  of  muscle 
fatigue  [18],  [14],  [12],  [25],  [29],  [11],  [13],  [27].  For  instance  it  may  be  possible  to  assess 
muscle  fatigue  by  determining  the  wavelet  decomposition  of  the  signal  with  the  wavelets  Sym4  or 
Sym5  and  eight  or  nine  levels  of  iterations  in  the  decomposition  [18].  In  addition,  some  authors 
have  proposed  to  estimate  the  MNF  of  a  SEMG  signal  by  using  WP.  By  comparing  (1)  MNF 
curves  derived  from  WP  with  (2)  MNF  curves  obtained  from  STFT,  it  is  inferred  that  both  (1)  and 
(2)  exhibit  qualitatively  similar  shapes  and  trends  [25].  WP  does  not  show  cross  terms  as  does  the 
WVD  and  has  the  advantage,  compared  to  the  traditional  TF  distributions,  of  being  the  only  one 
that  allows  selecting  the  frequency  range  of  interest  [29]. 

A  comparison  of  the  STFT,  the  CWD,  the  WVD,  and  the  continuous  wavelet  transform  (CWT)  is 
performed  in  [11]  and  [13].  In  this  comparison  it  is  inferred  that  the  STFT  is  affected  in  time  or  in 
frequency  resolution  by  the  choice  of  the  time  window.  It  is  also  concluded  that  the  WVD  is 
disturbed  by  cross  terms  affecting  the  estimation  of  the  frequency.  Moreover,  the  CWD  is  the 
most  suitable  method  to  represent  non-stationary  SEMG  signals;  however,  its  use  requires  careful 
selection  of  the  parameters  in  the  kernel  which  should  be  adapted  to  each  specific  signal.  Finally, 
the  CWT  is  very  reliable  for  the  analysis  of  non-stationary  signals,  the  IF  estimates  based  on  the 
CWT  are  less  noisy  than  the  other  techniques,  and  the  CWT  does  not  require  any  smoothing 
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function.  The  main  drawback  of  the  CWT  is  that  it  is  computationally  expensive  and  requires 
large  storage  requirements.  Discrete  wavelet  transform  (DWT)  is  proposed  to  overcome  this  issue 

[18] ,  [14],  [12],  [25],  [29], 

All  the  analyzed  papers  regarding  the  assessment  of  muscle  fatigue  by  using  wavelets  and  WP 
have  focused  in  analyzing  just  changes  in  the  frequency  of  the  SEMG  signal  over  time,  yet  the 
amplitude  of  the  SEMG  signal  has  been  disregarded. 

One  of  the  characteristics  of  muscle  fatigue  is  a  decrease  in  the  frequency  and  an  increase  in  the 
amplitude  of  the  SEMG  signal.  By  using  a  linear  regression  model,  slopes  can  be  computed  to 
measure  the  decrease  or  increase  rate  of  the  estimated  frequency  and  amplitude  over  time  [15], 
[20],  [30],  [19],  [27],  [32],  [9].  A  joint  analysis  of  frequency  and  amplitude  is  proposed  to  classify 
the  slopes  into  four  stages:  (1)  increase  force  when  both  the  frequency  and  the  amplitude  slopes 
increase  over  time,  (2)  recovery  when  the  frequency  slopes  increase  and  the  amplitude  slopes 
decrease,  (3)  decrease  force  when  both  the  frequency  and  the  amplitude  slopes  decrease  over 
time,  or  (4)  fatigue  when  the  frequency  slopes  decrease  and  the  amplitude  slopes  increase  [20], 

[19] .  Although  a  joint  analysis  of  frequency  and  amplitude  seems  to  be  suitable  to  identify  muscle 
fatigue,  a  relationship  is  not  yet  known  between  the  decrease  rate  of  the  frequency  and  the 
increase  rate  of  the  amplitude,  and  the  intensity  of  fatigue  that  is  felt  by  the  subject.  Additionally, 
the  time  interval  to  compute  the  slopes  is  also  unknown.  Some  papers  compute  the  slopes  from 
the  whole  set  of  the  estimated  IF  and  IA  while  others  compute  the  slopes  on  segments  of  20s,  40s 
or  50s  of  the  estimated  IF  and  IA  [30],  [19],  [27],  [32],  [9]. 

1.3  Desired  Characteristics  of  a  Time-Frequency  (TF) 
Distribution  in  the  Analysis  of  SEMG  Signals 

Time-Frequency  analysis  is  an  important  tool  in  the  development  of  indices  of  muscle  fatigue 
since  the  frequency  and  the  amplitude  of  the  SEMG  signal  can  be  extracted  from  a  TF 
distribution.  The  analyzed  papers  have  focused  on  the  estimation  of  the  frequency,  yet  they  have 
disregarded  the  evolving  of  the  amplitude  over  time.  In  this  research  both  the  frequency  and 
amplitude  are  extracted  from  the  spectrogram,  the  CWD,  and  the  SPWVD.  The  evaluation  of 
these  techniques  is  based  on  the  desired  characteristics  of  a  TF  distribution  mentioned  below. 

1.3.1  Time  and  Frequency  Marginals 

The  time  and  frequency  marginals  of  a  joint  TF  distribution  P(t,  u)  of  a  signal  x(t)  are 
respectively  given  by 

/OO 

P(t,Lu)df  =  \x(t)\2  (1.5) 

-OO 
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P(t,u)dt  =  |X(w)|2 

-OO 


x(t)e  juJtdt 


(1.6) 


where  X(u)  is  the  Fourier  transform  of  the  signal  x(t),  P(t)  is  the  time  marginal,  and  P( u)  is  the 
frequency  marginal  of  the  signal  x(t).  If  the  TF  distribution  satisfies  the  time  and  frequency 
marginals,  then  according  to  (1.1)  and  (1.2)  the  time  and  frequency  marginals  must  be, 
respectively,  the  instantaneous  energy  and  the  energy  density  spectrum  of  the  analyzed  signal. 


1.3.2  Time  and  Frequency  Finite  Support 

If  a  signal  x(t)  is  zero  inside  the  time  interval  [t\,  t2\,  then  a  TF  distribution  satisfying  P(t,  u)  =  0 
for  t  G  [t],  t2]  is  expected.  Similarly,  if  the  spectrum  of  a  signal  x(  t)  is  zero  inside  the  frequency 
interval  [oq.o^],  then  aTF  distribution  satisfying  P(t,u)  =  0  force  G  [u>i.u>2]  is  expected  [6],  [23]. 

1.3.3  Time  and  Frequency  Localization 

The  time  localization  of  a  signal  x(t)  at  time  t0  is  given  by 

x(t)  =  5(t  —  to)  =>•  P(t,u)  =  5(t  —  t0), 

and  the  frequency  localization  at  frequency  uo  is  given  by 

X{u)  =  5{u  —  ceo)  P(t,  u)  =  5(co  —  ceo) 


[23], 

1.3.4  Positivity 

Wigner  demonstrated  that  any  bilinear  distribution  that  satisfies  the  time  and  frequency  marginals 
cannot  always  be  positive;  it  must  have  regions  of  negative  values  for  any  signal.  However,  Cohen 
showed  that  there  exist  TF  distributions  that  are  not  bilinear  but  are  positive  and  satisfy  the 
marginals.  One  of  the  advantages  of  having  a  TF  distribution  that  is  positive  and  satisfies  the 
marginals  is  that  the  distribution  also  satisfies  the  time  and  frequency  finite  support  [6]. 

1.3.5  Suppression  of  Cross  Terms 

Cross  terms  are  artifacts  or  interference  terms  that  are  undesirable  in  the  analysis  of  muscle 
fatigue  because  they  negatively  affect  the  estimation  of  the  IF.  These  cross  terms  are  introduced 
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by  the  bilinearity  of  some  distributions  such  as  the  Wigner-Ville  distribution.  For  instance,  let  a 
signal  x(t)  be  expressed  as  the  sum  of  two  parts  x\ (t)  +  x2(t).  The  WVD  of  x(t)  is  then  given  by 

WVD(t,u)  =  WVDu(t,w)  +  WVD22(t,u)  +  WVD12(t,  u)  +  WVD21(t,u), 

where  WV  D\2(t,u)  and  WV  1)2\ (t.uo)  are  the  cross-WVD.  Therefore,  the  WVD  of  the  sum  of 
two  signals  is  not  the  sum  of  the  WVD  of  each  signal.  It  also  adds  cross  terms  interfering  with  the 
analysis  of  the  distribution. 


1.3.6  Window  Independence 

The  aim  of  introducing  a  TF  window  is  to  reduce  cross  terms  and  better  localize  the  quick  changes 
of  the  SEMG  signal.  However,  these  windows  destroy  desirable  properties  of  a  distribution  such 
as  the  time  and  frequency  marginals  and  time  and  frequency  support.  Furthermore,  the  choice  of 
the  size  of  the  TF  windows  play  an  important  role  in  the  assumption  of  the  stationarity  of  the 
SEMG  signal  as  well  as  in  a  TF  localization  trade-off.  A  TF  localization  trade-off  implies  that  a 
good  localization  in  time  leads  to  a  poor  localization  in  frequency  and  vice-versa. 


1.3.7  Energy  Conservation 

The  total  energy  of  an  ideal  joint  TF  distribution  should  be  the  total  energy  of  the  signal  x(t)  [6] 

/oo  r  oo  r  oo  r  oo 

/  P(t,u>)dwdt  =  /  \x(t)\2dt  =  /  \X(co)\2du. 

-oo  J — OO  J  —  oo  J — oo 


1.3.8  Instantaneous  Frequency 

The  first  order  moment  of  frequency  for  a  given  time  defines  the  IF  of  a  signal  x(t)  as 

fZoP(t,Uj)duJ  ' 

If  the  signal  is  written  in  terms  of  its  phase  <p  and  amplitude  A, 

x(t)  =  A(t)e?vW, 


(1.7) 


then  it  is  straightforward  to  show  that 


4>{t) 


Id 
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(1.8) 
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1.4  TF  distributions 


The  following  TF  distributions  are  the  most  widely  used  in  the  analysis  of  muscle  fatigue. 

1.4.1  The  Short-Time  Fourier  Transform  (STFT)  and  the  Spectrogram 

The  linear  STFT  is  defined  as  the  Fourier  transform  of  subsequences  (short  segments)  of  a  signal 
windowed  in  the  time  domain.  As  a  result  the  STFT  is  a  function  of  the  time  t  and  the  frequency 
u  depending  of  the  size  of  the  window  h  (t  ). 

STFT(t,  u)  =  J  x(r)h(T  -  t)e~j2™TdT, 

where  t  is  the  shift  of  the  time  window.  A  time  Gaussian  window  of  size  128  samples  is  used  in 
this  research. 

The  spectrogram  (squared  magnitude  of  the  STFT)  is  always  positive,  real,  and  computationally 
fast.  However,  due  to  the  time  window  introduced  by  the  STFT,  the  spectrogram  does  not  satisfy 
the  marginals  as  well  as  the  TF  support  and  the  TF  localization.  The  spectrogram  suffers  from  a 
TF  localization  trade-off  and  assumes  stationarities  of  the  signal  in  short  intervals  of  time  given 
by  the  size  of  the  window  h(t).  The  previous  assumption  of  stationarity  is  true  for  SEMG  signals 
with  static  contractions. 

1.4.2  The  Wigner-Ville  Distribution  (WVD) 

The  bilinear  WVD  distribution  is  defined  as: 

W(t,u)  =  ±.Jx-(t-1-T)x(t  +  It)  e-^dr, 

where  *  denotes  the  complex  conjugate  of  a  function.  This  distribution  differs  from  the 
spectrogram  in  that:  (1)  it  is  not  always  positive,  (2)  it  holds  the  time  and  frequency  marginals,  (3) 
it  generates  considerable  amount  of  cross  terms  in  the  analysis  of  multicomponent  signals,  and  (4) 
it  is  non  local.  The  WVD,  like  the  spectrogram,  does  not  satisfy  the  time  and  frequency  finite 
support.  In  addition,  the  WVD  has  not  been  used  in  the  analysis  of  muscle  fatigue  since  cross 
terms  introduce  errors  in  the  estimation  of  the  IF. 

A  bidimensional  smoothing  function  (kernel)  can  be  designed  to  convolve  with  the  WVD.  The 
aim  of  introducing  this  kernel  is  to  make  the  WVD  more  local  and  to  reduce  the  cross  terms. 
Additionally,  desirable  properties  can  be  accomplished  by  constraining  the  kernel  function.  The 
CWD  and  the  SPWVD  define  two  different  kernels  to  make  the  WVD  suitable  in  the  analysis  of 
muscle  fatigue. 
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1.4.3  The  Choi- Williams  distribution  (CWD) 

The  kernel  for  the  CWD  is  defined  as 


(7 tut) 

=  e~  2^  .  (1.9) 

It  is  noticed  that  for  this  distribution  when  a  — »  +00,  the  WVD  is  obtained.  Likewise,  the  cross 
terms  are  almost  zero  when  a  — >  0.  The  CWD  satisfies  the  time  and  frequency  marginals  and 
localization,  but  it  is  still  non  positive  and  it  does  not  satisfy  the  time  and  frequency  finite  support. 
However,  its  efficiency  strongly  depends  on  the  choice  of  a  critical  parameter  which  has  to  be 
adapted  according  to  the  nature  of  the  SEMG  signal.  The  parameter  a  was  set  to  1,  and  time  and 
frequency  Gaussian  windows  were  set  to  128  and  64  samples,  respectively,  for  this  research. 

1.4.4  The  Smoothed  Pseudo  Wigner-Ville  Distribution  (SPWVD) 

It  was  shown  before  that  the  spectrogram  is  dependent  on  a  short-time  window  yielding  to  a  TF 
localization  trade-off.  In  order  to  make  the  smoothing  function  more  independent  in  both  time 
and  frequency,  the  SPWVD  introduces  a  separable  kernel  function  defined  as 

=  g(r)H(u), 

where  H{y)  is  a  frequency  smoothing  window  with  H( 0)  being  forced  to  1,  and  g{r)  is  a 
smoothing  time  window  with  <7(0)  being  forced  to  1.  However,  because  of  the  introduction  of  the 
kernel,  the  time  and  frequency  marginals  of  the  WVD  are  not  satisfied.  The  SPWVD  is  still  non 
positive  and  does  not  hold  the  time  and  frequency  finite  support.  Time  and  frequency  Gaussian 
windows  of  lengths  128  and  64  samples,  respectively,  are  used  for  this  research. 

1.5  Evaluation  of  the  Time-Frequency  (TF)  Distributions 

In  this  section  a  comparison  of  the  spectrogram,  WVD,  CWD  and  SPWVD  is  made  by  estimating 
known  parameters,  the  IF,  and  the  IA  of  synthesized  signals.  Matlab  functions  from  the  TF 
toolbox  downloaded  from  tftb.nongnu.org.  are  used  to  compute  each  TF  distribution. 

1.5.1  Tests  with  Known  signals 

The  following  eight  experiments  use  known  signals  whose  IF  and  IA  values  can  be  theoretically 
computed.  These  theoretical  IF  and  IA  values  are  compared  with  those  values  estimated  from  the 
Spectrogram,  the  WVD,  the  CWD,  and  the  SPWVD.  Based  on  these  comparisons,  the  techniques 
that  better  estimate  the  IF  and  IA  of  known  signals  are  selected  to  estimate  the  IF  and  IA  of 
SEMG  signals.  The  results  of  each  experiment  are  shown  at  the  end  of  this  Chapter. 
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Experiment  1 

Two  sinusoids  of  different  frequencies  modulated  in  amplitude 

Consider  a  signal  which  is  the  sum  of  two  sinusoids  of  different  frequencies  and  is  modulated  in 
amplitude  by  multiplying  it  by  a  ramp  signal. 

s(t)  =  At(sm(2n  fit  +  sin(27r/2f))  (1-10) 

where  A  is  equal  to  1,  f\  =  100  Hz  and  /2  =  150  Hz.  Knowing  the  frequency  components  of  s(t) 
is  easy  to  compute  the  IF  as  /  =  =  125  Hz.  The  expected  IA  is  a  ramp  function.  See  Figs. 

1.1  and  1.2.  The  Matlab  code  to  generate  this  signal  is  shown  below. 

Fs=1024;  %sampling  frequency 
fl=100;  %  frequency  of  the  first  sinusoid 
f2=150;  %  frequency  of  the  second  sinusoid 
n=0:l/ (Fs-1) :1;  %  number  of  samples 
sl=sin (2*pi*f l*n) ;  %  first  sinosoid 
s2=sin (2*pi*f2*n) ;  %  second  sinusoid 

s3=(sl+s2) .  *n;  %  sum  of  two  sinusoids  modulated  in  amplitude 


Experiment  2 

Two  parallel  chirp  signals  increasing  in  frequency  and  modulated  in  amplitude 

A  chirp  signal  is  characterized  by  its  frequency  content  changing  with  time  and  therefore,  it  is  a 
non- stationary  signal  which  is  very  important  to  analyze.  Using  the  command  chirp  in  Matlab, 
a  composed  signal  is  generated  by  adding  two  chirp  signals.  The  first  chirp  signal  has  an  IF  going 
from  100  Hz  to  150  Hz,  and  the  second  chirp  signal  has  an  IF  going  from  125  Hz  to  175  Hz,  both 
with  1024  samples.  The  expected  instantaneous  frequency  of  the  added  signals  starts  at  112.5  Hz 
and  ends  at  162.5  Hz.  This  signal  is  modulated  in  amplitude  by  multiplying  it  by  a  Gaussian 
signal.  See  Figs.  1.3  and  1.4.  The  Matlab  code  to  generate  this  signal  is  shown  below. 

Fs=1024;  %  sampling  frequency 
n=0:l/  (Fs-1)  :1;  %  number  of  samples 

amp=amgauss (1024, 512, 256) ;  %  Gaussian  signal  centralized  at 
time  sample  512 

sl=chirp (n, 100, 1, 150) ;  %  first  chirp  signal 
s2=chirp  (n,  125,  1,  175) ;  %  second  chirp  signal 

s=(sl+s2) ,*amp;  %  sum  of  two  chirp  signals  modulated  in  amplitude 


Experiment  3 

Two  crossed  chirp  signals  increasing  and  decreasing  in  frequency,  and  modulated  in  amplitude 

Using  the  command  chirp  in  Matlab,  a  composed  signal  is  generated  adding  two  chirp  signals. 
The  first  chirp  signal  has  an  IF  going  from  100  Hz  to  300  Hz  and  the  second  chirp  signal  has  an  IF 
going  from  300  Hz  to  100  Hz,  both  with  1024  samples.  The  expected  IF  of  the  composed  signal  is 
a  constant  with  a  value  of  200  Hz.  This  signal  is  modulated  in  amplitude  by  multiplying  it  by  a 
Gaussian  signal.  See  Figs.  1.5  and  1.6.  The  Matlab  code  to  generate  this  signal  is  shown  below. 
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Fs=1024;  %  sampling  frequency 
n=0:l/  (Fs-1)  :1;  %  number  of  samples 

amp=amgauss (1024, 512, 256) ;  %  Gaussian  signal  centralized  at 
time  sample  512 

sl=chirp (n, 100, 1, 300) ;  %  first  chirp  signal 
s2=chirp (n, 300, 1, 100) ;  %  second  chirp  signal 

s=(sl+s2) . *amp;  %  sum  of  two  chirp  signals  modulated  in  amplitude 


Experiment  4 

Chirp  signal  with  increasing  frequencies  followed  by  a  sinusoid  and  a  Chirp  signal  with 
decreasing  frequencies 

Another  chirp  signal  is  computed  with  the  following  characteristics:  (1)  constant  amplitude  and 
linear  frequency  modulation  varying  from  100  Hz  to  400  Hz  in  the  interval  [0,  384]  samples,  (2) 
constant  amplitude  and  constant  frequency  of  400  Hz  in  the  interval  [385,  640]  samples,  and  (3) 
constant  amplitude  and  linear  frequency  modulation  varying  from  400  Hz  to  100  Hz  in  the 
interval  [641, 1024]  samples.  See  Figs.  1.7  and  1.8.  The  Matlab  code  to  generate  this  signal  is 
shown  below. 

Fs=1024;  %  sampling  frequency 
n=0:l/  (Fs-1)  :1;  %  number  of  samples 

sl=fmlin ( 384 , 0 . 1 , 0 . 4 ) ;  %  first  signal  linearly  modulated  in  frequency 
s2=fmlin (256, 0 . 4, 0 . 4 ) ;  %  second  signal  linearly  modulated  in  frequency 
s3=fmlin (384, 0 . 4, 0 . 1) ;  %  third  signal  linearly  modulated  in  frequency 
s= [ si ; s2 ; s3 ] ;  %  three  signals  modulated  in  frequency 


Experiment  5 

Four  successive  chirp  signals  increasing  in  frequency 

Consider  a  signal  with  constant  amplitude  and  linear  frequency  modulation  varying  from  100  Hz 
to  400  Hz  in  four  time  intervals  of  256  samples  each  one.  See  Figs.  1.9  and  1.10.  The  Matlab 
code  to  generate  this  signal  is  shown  below. 

Fs=1024;  %  sampling  frequency 
n=0:l/  (Fs-1)  :1;  %  number  of  samples 

sl=fmlin (256, 0 . 1, 0 . 4) ;  %  first  signal  linearly  modulated  in  frequency 
s2=fmlin (256, 0 . 1, 0 . 4) ;  %  second  signal  linearly  modulated  in  frequency 
s3=fmlin (256, 0 . 1, 0 . 4) ;  %  third  signal  linearly  modulated  in  frequency 
s4=fmlin (256, 0 . 1, 0 . 4) ;  %  fourth  signal  linearly  modulated  in  frequency 
s= [ si ; s2 ; s3 ; s4 ] ;  %  four  signals  modulated  in  frequency 
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Experiment  6 

Three  atoms  localized  in  frequency  and  time 

An  atom  is  an  elementary  waveform  localized  at  a  certain  time  and  frequency  over  the  TF  plane. 
A  linear  combination  of  elementary  atoms  using  the  Matlab  function  atoms  from  the 
time-frequency  toolbox  downloaded  from  tftb.nongnu.org  is  generated.  This  experiment  uses 
three  atoms  centralized  in  time  and  frequency  according  to:  (1)  atom  1  is  centralized  in  time  in 
the  sample  64  with  a  duration  of  100  samples  and  in  frequency  200  Hz,  (2)  atom  2  is  centralized 
in  time  in  the  sample  320  with  a  duration  of  100  samples  and  in  frequency  400  Hz,  and  finally  (3) 
atom  3  is  centralized  in  time  in  the  sample  640  with  a  duration  of  100  samples  and  in  frequency 
250  Hz.  See  Table  1.1  and  Figs.  1.12  and  1.13.  The  ideal  TF  distribution  of  these  three  atoms  is 
represented  in  the  TF  plane  shown  in  Fig.  1.11. 

Table  1.1  Description  of  the  time-frequency  localization  of  three  atoms  used  in  Experiment  6. 


Coordinates 

Atoms 

Time  (samples) 

Frequency  (Hz) 

Duration  (samples) 

1 

64 

200 

100 

2 

320 

400 

100 

3 

640 

250 

100 

The  Matlab  code  to  generate  this  signal  is  shown  below. 


s=atoms (1024, 


[  64, 0. 2, 100, 1;  32 0, 0. 4, 100, 1;  64 0,0. 25, 100,1]); 


Experiment  7 

Four  atoms  localized  in  frequency  and  time 

This  experiment  uses  four  atoms  centralized  in  time  and  frequency  according  to:  (1)  atom  1  is 
centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in  frequency  1 10  Hz,  (2) 
atom  2  is  centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in  frequency 
330  Hz,  (3)  atom  3  is  centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in 
frequency  420  Hz,  and  finally  (4)  atom  4  is  centralized  in  time  in  the  sample  512  with  a  duration 
of  200  samples  and  in  frequency  200  Hz.  See  Table  1.2  and  Figs.  1.15  and  1.16.  The  ideal  TF 
distribution  of  these  four  atoms  is  represented  in  the  TF  plane  shown  in  Fig.  1.14. 

The  Matlab  code  to  generate  this  signal  is  shown  below. 

s=atoms (1024, [200, 0.11, 100,1; 200, 0.33, 100, 1; 200, 0.42, 100,1 ; 

512, 0.2,200, 1] ) ; 
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Table  1.2  Description  of  the  time-frequency  localization  of  four  atoms  used  in  Experiment  7. 


Coordinates 

Atoms 

Time  (samples) 

Frequency  (Hz) 

Duration  (samples) 

1 

200 

110 

100 

2 

200 

330 

100 

3 

200 

420 

100 

4 

512 

200 

200 

Experiment  8 

Five  atoms  localized  in  frequency  and  time 

This  experiment  uses  five  atoms  centralized  in  time  and  frequency  according  to:  (1)  atom  1  is 
centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in  frequency  110  Hz,  (2) 
atom  2  is  centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in  frequency 
330  Hz,  (3)  atom  3  is  centralized  in  time  in  the  sample  200  with  a  duration  of  100  samples  and  in 
frequency  420  Hz,  (4)  atom  4  is  centralized  in  time  in  the  sample  512  with  a  duration  of  200 
samples  and  in  frequency  200  Hz,  and  finally  (5)  atom  5  is  centralized  in  time  in  the  sample  615 
with  a  duration  of  25  samples  and  in  frequency  450  Hz.  See  Table  1.3  and  Figs.  1.18  and  1.19. 
The  ideal  TF  distribution  of  theses  five  atoms  is  represented  in  the  TF  plane  shown  in  Fig.  1.17. 


Table  1.3  Description  of  the  time-frequency  localization  of  five  atoms  used  in  Experiment  8. 


Coordinates 

Atoms 

Time  (samples) 

Frequency  (Hz) 

Duration  (samples) 

1 

200 

110 

100 

2 

200 

330 

100 

3 

200 

420 

100 

4 

512 

200 

200 

5 

615 

450 

25 

The  Matlab  code  to  generate  this  signal  is  shown  below. 


s=atoms (1024, [200, 0.11, 100,1; 200, 0.33, 100,1; 200, 0.42, 100,1 ; 
512, 0.2,200, 1] ) ; 
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1.6  Evaluation  of  the  Time-Frequency  (TF)  Distributions 


This  section  presents  the  results  after  having  compared  the  IF  and  IA  values  computed 
theoretically  from  the  known  signals  with  those  IF  and  IA  values  estimated  from  the 
Spectrogram,  the  WVD,  the  CWD,  and  the  SPWVD.  According  to  Figs.  1.1,  1.3,  1.5,  1.7,  1.9, 
1.12,  1.15,  and  1.18  the  WVD  is  not  window  dependent  and  gives  the  best  TF  resolution,  but 
strongly  presents  cross  terms.  The  spectrogram  is  time-window  dependent  and  gives  the  worst  TF 
resolutions  but  with  negligible  cross  terms.  Also,  the  SPWVD  is  window  independent  but  allows 
to  choose  the  best  compromise  between  resolution  and  cross  terms.  In  the  theoretical 
experiments,  it  is  shown  in  Figs.  1.2,  1.4,  1.6,  1.8,  1.10,  1.13,  1.16,  and  1.19  that  although  the 
spectrogram  and  the  SPWVD  are  window  dependent  and  do  not  satisfy  the  TF  finite  support,  they 
can  give  good  estimations  of  the  IF  and  the  IA.  In  contrast  the  WVD  and  the  CWD  provide 
negative  frequencies  and  assume  frequencies  that  are  outside  the  bandwidth  of  the  analyzed 
signal.  Table  1.4  summarizes  the  strengths  and  shortcomings  of  each  evaluated  TF  distribution 
(spectrogram,  WVD,  CWD,  SPWVD). 


Table  1.4  Summary  of  strengths  and  shortcomings  of  the  evaluated  techniques. 


Spectrogram 

WVD 

CWD 

SPWVD 

Advantages 

Fast  to  compute 

Energy  conservation 

Energy  conservation 

Reduces  cross  terms 

Positive 

TF  marginals 

TF  marginals 

More  local 

TF  localization 

TF  localization 

Finite  time  support 

Reduces  cross  terms 

Non-local 

More  local 

Disadvantages 


TF  localization  trade-off 

Non-positive 

Non-positive 

Non-positive 

Not  satisfy  TF  marginals 

Cross  terms 

Not  satisfy  TF  support 

TF  localization  trade-off 

Not  satisfy  TF  support 

Negative  IF 

Negative  IF 

Not  satisfy  TF  marginals 

Assumes  signal  stationary 

Not  satisfy  TF  support 
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Figure  1.1  TF  distributions  of  the  signal  in  Experiment  1 .  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  a  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.2  Computed  IF  and  IA  for  Experiment  1.  The  expected  IF  is  125  Hz,  and  the  expected 
IA  is  a  ramp  function.  The  WVD  and  the  CWD  provide  negative  IF  and  frequencies  outside  the 
bandwidth  of  the  signal,  (a)  The  spectrogram  gives  a  good  estimation  of  the  IF  and  the  IA.  (b) 
The  WVD  gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the  IA.  (c)  The  CWD 
gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the  IA.  (d)  The  SPWVD  gives  a  good 
estimation  of  the  IF  and  the  IA. 
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Figure  1.3  TF  distributions  of  the  signal  in  Experiment  2.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  The  region  of  support  of  the  composed  signal  is  in  the  time  sample 
[384,640].  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b)  WVD.  It 
gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian  window 
of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter  a  =  1,  as 
defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and  frequency 
Gaussian  window  of  length  64  samples). 
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Figure  1.4  Computed  IF  and  IA  for  Experiment  2.  The  expected  IF  is  1 12.5  Hz  at  time  sample 
1  and  162.5  Hz  at  time  sample  1024,  and  the  expected  IA  is  a  Gaussian  function  centered  at 
time  sample  512  with  duration  256  samples.  The  WVD  and  the  CWD  provide  negative  IF  and 
frequencies  outside  the  bandwidth  of  the  signal.  Any  of  the  evaluated  TF  distributions  satisfy  the 
TF  support  since  all  of  them  show  IF  outside  the  region  of  support,  (a)  The  spectrogram  gives  a 
good  estimation  of  the  IF  and  the  IA.  (b)  The  WVD  gives  a  wrong  estimation  of  the  IF  and  a  good 
estimation  of  the  IA.  (c)  The  CWD  gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the 
IA.  (d)  The  SPWVD  gives  a  good  estimation  of  the  IF  and  the  IA. 
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Figure  1.5  TF  distributions  of  the  signal  in  Experiment  3.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  The  region  of  support  of  the  composed  signal  is  in  the  time  sample 
[384,640].  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b)  WVD.  It 
gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian  window 
of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter  a  =  1,  as 
defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and  frequency 
Gaussian  window  of  length  64  samples. 
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Figure  1.6  Computed  IF  and  IA  for  Experiment  3.  The  expected  IF  is  a  constant  frequency  of 
200  Hz  along  the  time,  and  the  expected  IA  is  a  Gaussian  function  centered  at  time  sample  512 
with  duration  256  samples.  The  WVD  and  the  CWD  provide  negative  IF  and  frequencies  outside 
the  bandwidth  of  the  signal.  Any  of  the  evaluated  TF  distributions  satisfy  the  TF  support  since  all 
of  them  show  IF  outside  the  region  of  support,  (a)  The  spectrogram  gives  a  good  estimation  of  the 
IF  and  the  IA.  (b)  The  WVD  gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the  IA. 
(c)  The  CWD  gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the  IA.  (d)  The  SPWVD 
gives  a  good  estimation  of  the  IF  and  the  IA. 
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Figure  1.7  TF  distributions  of  the  signal  in  Experiment  4.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.8  Computed  IF  and  IA  for  Experiment  4.  First,  the  expected  IF  increases  linearly  from 
100  Hz  to  400  Hz  at  time  sample  [0, 384].  Second,  the  expected  IF  is  a  constant  of  400  Hz  at  time 
sample  [385,  640] .  Finally,  the  expected  IF  decreases  linearly  from  400  Hz  to  100  Hz  at  time  sample 
[641, 1024],  The  expected  IA  is  a  constant  value.  The  WVD  and  the  CWD  provide  frequencies 
outside  the  bandwidth  of  the  signal.  All  the  evaluated  TF  distributions  show  transitions  in  the 
discontinuities  of  the  analyzed  signal  giving  a  wrong  estimation  of  the  IA.  These  abrupt  changes 
of  the  signal  are  not  presented  in  SEMG  signals,  (a)  The  spectrogram  gives  a  good  estimation  of 
the  IF.  (b)  The  WVD  gives  a  wrong  estimation  of  the  IF.  (c)  The  CWD  gives  a  wrong  estimation 
of  the  IF.  (d)  The  SPWVD  gives  a  good  estimation  of  the  IF. 
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Figure  1.9  TF  distributions  of  the  signal  in  Experiment  5.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.10  Computed  IF  and  IA  for  Experiment  5.  The  expected  IF  increases  linearly  from 
100  Hz  to  400  Hz  in  four  time  intervals  of  256  samples.  The  expected  IA  is  a  constant  value. 
The  spectrogram  and  the  SPWVD  are  less  susceptible  to  the  discontinuities  of  the  signal  giving  a 
smooth  IA.  (a)  The  computed  IF  from  the  spectrogram  does  not  reach  the  maximum  frequency  of 
400  Hz.  (b)  The  WVD  gives  a  good  estimation  of  the  IF.  (c)  The  CWD  gives  a  good  estimation  of 
the  IF.  (d)  The  computed  IF  from  SPWVD  also  does  not  reach  the  maximum  frequency  of  400  Hz. 
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Figure  1.11  Ideal  TF  distribution  of  the  signal  in  Experiment  6  with  a  normalized  frequency  from 
0  to  0.5. 
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Signal 


Figure  1.12  TF  distributions  of  the  signal  in  Experiment  6.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.13  Computed  IF  and  IA  for  Experiment  6.  Three  Gaussian  atoms  in  the  TF  plane.  One 
atom  is  located  at  time  sample  64  with  an  expected  IF  of  200  Hz,  a  second  atom  is  located  at  time 
sample  320  with  an  expected  IF  of  400  Hz,  and  a  third  atom  is  located  in  sample  640  with  an 
IF  of  250  Hz.  The  expected  IA  for  the  three  atoms  is  three  Gaussian  functions  centered  at  time 
samples  64,  320  and  640.  (a)  The  spectrogram  gives  a  good  estimation  of  the  IF  and  the  IA.  The 
spectrogram  is  the  only  one  that  shows  an  IF  of  zero  when  the  signal  is  zero.  In  this  case  the 
time  support  property  is  satisfied,  (b)  The  WVD  gives  a  wrong  estimation  of  the  IF  and  a  good 
estimation  of  the  IA.  (c)  The  CWD  gives  a  wrong  estimation  of  the  IF  and  a  good  estimation  of  the 
IA.  (d)  The  SPWVD  gives  a  good  estimation  of  the  IF  and  the  IA. 
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Figure  1.14  Ideal  TF  distribution  of  the  signal  in  Experiment  7  with  a  normalized  frequency  from 
0  to  0.5. 


28 


Signal 


Signal 


Signal 


5000  0  200  400  600  800  1000 

Time  Sample 


(c) 


(d) 


Figure  1.15  TF  distributions  of  the  signal  in  Experiment  7.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.16  Computed  IF  and  IA  for  Experiment  7.  Four  Gaussian  atoms  in  the  TF  plane.  The 
expected  IF  is  287  Hz  at  time  sample  [100, 300]  and  200  Hz  at  time  sample  [312,  711].  The  expected 
IA  for  the  four  atoms  is  two  Gaussian  functions  centered  at  time  sample  200  and  512.  (a)  The 
spectrogram  gives  a  good  estimation  of  the  IF  and  the  IA.  (b)  The  WVD  gives  a  wrong  estimation 
of  the  IF  and  a  good  estimation  of  the  IA.  (c)  The  CWD  gives  a  wrong  estimation  of  the  IF  and  a 
good  estimation  of  the  IA.  (d)  The  SPWVD  gives  a  good  estimation  of  the  IF  and  the  IA. 
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Figure  1.17  Ideal  TF  distribution  of  the  signal  in  Experiment  8  with  a  normalized  frequency  from 
0  to  0.5. 
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Figure  1.18  TF  distributions  of  the  signal  in  Experiment  8.  Each  TF  distribution  has  a  normalized 
frequency  from  0  to  0.5.  (a)  Spectrogram  with  a  time  Gaussian  window  of  length  128  samples,  (b) 
WVD.  It  gives  the  best  TF  resolution  but  is  affected  by  cross  terms,  (c)  CWD  with  a  time  Gaussian 
window  of  length  128  samples,  frequency  Gaussian  window  of  length  64  samples,  and  parameter 
a  —  1,  as  defined  in  (1.9).  (d)  SPWVD  with  a  time  Gaussian  window  of  length  128  samples  and 
frequency  Gaussian  window  of  length  64  samples. 
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Figure  1.19  Computed  IF  and  IA  for  Experiment  8.  Five  Gaussian  atoms  in  the  TF  plane.  The 
expected  IF  is  287  Hz  at  time  sample  [100,  300],  200  Hz  at  time  sample  [312,  711],  and  450  Hz  at 
time  sample  [590,  639].  The  expected  IA  for  the  five  atoms  is  three  Gaussian  functions  centered 
at  time  sample  200,  512,  and  615.  In  general  any  of  the  analyzed  TF  distributions  could  not 
reach  the  maximum  frequency  of  450  Hz.  The  WVD  and  the  CWD  still  provide  negative  IF 
and  frequencies  outside  the  bandwidth  of  the  signal.  The  spectrogram  and  the  SPWVD  provide 
adequate  estimations  of  the  IF  and  IA.  (a)  Spectrogram,  (b)  WVD.  (c)  CWD.  (d)  SPWVD. 
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2 


Materials  and  Method 

2.1  Materials 

Muscle  fatigue  was  evaluated  from  SEMG  signals  recorded  from  26  normal  human  subjects,  1 1 
females  and  15  males.  The  females  had  an  average  age  of  23  years,  weight  of  146  pounds  and 
height  of  163  cm;  the  males  had  an  average  age  of  27  years,  weight  of  184  pounds  and  height  of 
177  cm.  The  SEMG  signals  were  acquired  at  a  sampling  frequency  of  1024  Hz  from  six  neck 
muscles:  trapezius,  splenius  capitis  and  sternocleidomastoid,  from  both  left  and  right  sides. 

A  Delsys  Bagnoli-8  EMG  system  with  DE-2.1  standard  differential  EMG  electrodes  was  used  for 
the  collection  of  the  myoelectric  data.  Prior  to  entrance  into  the  study,  basic  anthropometry 
measurements  (age,  weight,  height)  were  collected  from  each  subject  to  ensure  that  they  fulfilled 
the  requirements  of  the  study. 

2.2  Study  Design 

For  this  research  3480  SEMG  recorded  signals  provided  by  the  Air  Force  Research  Laboratory 
are  analyzed  in  the  evaluation  of  muscle  fatigue.  The  objective  of  this  analysis  is  to  present 
indices  to  measure  neck  muscle  fatigue  during  prolonged  wear  of  weighted  flight  helmets.  These 
indices  are  useful  in  the  development  of  guidelines  to  guarantee  good  performance  and  safety  to 
the  pilots.  In  this  study  each  subject  was  exposed  to  five  sessions  of  eight  hours  with  a  minimum 
waiting  period  of  48  hours  between  them.  In  each  session  the  subject  wore  a  different  weighted 
flight  helmet  as  shown  in  Table  2.1.  They  also  performed  a  100%  maximum  voluntary  contraction 
(MVC)  before  and  after  the  eight-hour  session,  and  a  70%  MVC  at  the  end  of  every  hour  (first 
hour  to  seventh  hour).  The  70%  MVC  was  held  until  the  subject  could  no  longer  maintain  the 
exertion  for  a  maximum  time  of  three  minutes.  A  dynamometer  was  used  to  monitor  the  strength 
given  by  the  subject’s  neck  muscle  during  the  70%  and  100%  MVC  exertion.  During  this  strength 
monitoring,  SEMG  data  were  recorded  to  measure  the  change  in  muscular  activity  and  to  analyze 
muscle  fatigue.  During  each  session  of  eight  hours,  the  subjects  completed  a  comfort  survey 
before  the  initial  100%  MVC,  after  hours  two,  four  and  seven,  and  after  the  final  100%  MVC. 
After  collection  of  the  data,  the  SEMG  signal  of  each  subject  was  normalized  to  his  100%  MVC 
as  suggested  by  Merletti  [21]. 
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Table  2.1  Helmet  configurations  (Air  Force  Research  Laboratory). 


Cell 

Weight  (Lb) 

Helmet  Center  of  Gravity  (CG) 

A 

3.0 

Baseline 

B 

4.5 

Near  head 

C 

4.5 

Forward  head 

D 

6.0 

Near  head 

E 

6.0 

Forward  head 

2.3  Method 

The  idea  of  using  wavelets  in  the  analysis  of  muscle  fatigue  arises  as  an  alternative  solution  to  the 
TF  resolution  issue  (Heisenberg  uncertainty  principle)  introduced  by  the  STFT  and  SPWVD. 
Wavelets  suggest  that  this  resolution  problem  can  be  solved  by  defining  short  time  windows  for 
high  frequencies  and  long  time  windows  for  low  frequencies.  Additionally,  wavelets  have 
demonstrated  to  be  a  reliable  method  in  the  analysis  of  nonstationary  biological  signals  with  the 
advantage  of  being  independent  of  any  smoothing  function,  in  contrast  to  the  spectrogram, 
SPWVD,  and  CWD,  which  are  dependent  of  a  smoothing  function. 

The  proposed  technique  (see  Appendix  A)  to  assess  muscle  fatigue  starts  by  filtering  the  recorded 
SEMG  signal  with  a  cosine  modulated  filter  bank.  In  the  SEMG  tool  box,  a  cosine  modulated 
filter  bank  with  16  and  32  channels  is  provided.  The  outputs  of  the  16  or  32  channels  are  the 
subband  coefficients  which  are  then  squared  and  arranged  in  a  time  subband  representation.  The 
IF  was  estimated  from  Equation  1.7  as  the  first  order  moment  of  the  time  subband  representation 
at  a  certain  point  in  time.  To  estimate  the  IA  from  the  time  subband  representation,  the  following 
sequential  stages  are  proposed:  (1)  define  a  rectangular  window  for  each  subband  with  a  length 
that  is  inversely  proportional  to  the  central  frequency  of  the  subband;  (2)  find  in  each  subband  the 
maximum  coefficient  inside  the  window;  (3)  finally,  add  all  the  corresponding  maximum  values  to 
obtain  an  estimate  of  the  IA.  The  estimate  of  the  IA  for  the  next  point  in  time  is  computed  by 
shifting  the  rectangular  window  one  sample  to  the  right  and  then  steps  (1)  to  (3)  are  repeated  until 
the  IA  is  estimated  along  all  the  signal.  Table  2.2  summarizes  the  strengths  and  shortcomings  of 
using  a  filter  bank  as  a  TF  distribution  in  the  analysis  of  muscle  fatigue. 

After  the  IF  and  the  IA  are  estimated  from  the  signal,  IF  and  IA  slopes  are  computed  by  applying 
a  first  order  linear  regression  model  to  the  estimated  IF  and  IA  for  one  minute  intervals  every  one 
second.  The  aim  of  these  slopes  is  to  measure  the  decrease  or  increase  rate  of  the  IF  and  IA,  thus 
providing  an  index  to  measure  the  intensity  of  muscle  increase  force,  recovery,  muscle  decrease 
force  or  fatigue.  However,  it  was  not  found  in  the  literature  any  suggestion  about  the  decrease  or 
increase  rate  of  the  slopes  where  muscle  increase  force,  recovery,  muscle  decrease  force  and 
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Table  2.2  Strengths  and  shortcomings  of  using  a  filter  bank  as  a  TF  distribution. 


Advantages 

Positivity 

Good  IF  estimation 
Good  IA  estimation 
Reduced  cross  terms 
Window  independent 

Allow  to  select  the  frequency  range  of  interest 
Local  in  time  and  frequency 
Fast  to  compute 

Disadvantages 
TF  marginals  not  satisfied 
TF  support  not  satisfied 
No  energy  conservation 


fatigue  should  be  considered.  Likewise,  it  was  not  found  any  suggestion  about  either  the  time 
interval  to  be  used  in  the  first  order  linear  regression  model  or  the  time  shift  between  two 
consecutive  slopes  (see  Fig.  2.1). 

The  IF  and  the  IA  variables  estimated  from  the  filter  bank  are  compared  with  the  IF  and  the  IA 
variables  estimated  from  two  well  known  TF  distributions  used  in  the  analysis  of  biomedical 
signals,  namely,  spectrogram,  and  SPWVD.  First,  the  SEMG  signal  is  Hilbert  transformed  and 
then  any  of  the  mentioned  TF  distributions  is  computed.  Second,  the  IF  parameter  is  estimated  by 
computing  the  first  order  moment  of  the  frequency  distribution  at  a  certain  point  in  time  as  shown 
in  (1.7)  [15],  [21].  Finally,  the  IA  parameter  is  estimated  by  computing  the  square  root  of  the  time 
marginal  of  the  employed  TF  distribution.  To  compute  the  IF  and  the  IA  slopes,  a  first  order  linear 
regression  model  is  applied  to  the  estimated  IF  and  IA  with  the  same  time  intervals  and  time  shift 
between  slopes  used  in  the  filter  bank. 

The  IF  and  IA  slopes  derived  from  the  spectrogram  and  the  SPWVD  are  compared  and  correlated 
to  those  slopes  derived  from  the  filter  bank  by  using  a  relative  quadratic  difference  and  Pearson’s 
correlation  coefficient  [27].  The  IF  and  the  IA  of  the  SEMG  signal  are  not  estimated  by  using  the 
CWD  because  as  it  was  shown  in  Section  1.5  the  CWD  provides  an  IF  with  negative  frequencies 
and  frequencies  outside  the  bandwidth  of  the  signal. 

A  joint  analysis  of  frequency  and  amplitude  is  adopted  in  this  research  to  analyze  both  the 
decrease  or  increase  rate  of  the  IF  and  the  decrease  or  increase  rate  of  the  IA.  By  using  a  joint 
analysis  of  frequency  and  amplitude,  as  shown  in  Fig  2.2,  the  IF  and  IA  slopes  are  classified  into 
four  stages:  (1)  increase  force  when  both  the  IF  slopes  and  the  IA  slopes  increase  over  time,  (2) 
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Figure  2.1  The  time  interval  to  be  used  in  the  first  order  linear  regression  and  the  time  shift 
between  slopes  are  unknown  in  the  literature. 

recovery  from  previous  muscle  fatigue  when  there  is  an  increase  in  the  IF  slopes  and  a  decrease  in 
the  IA  slopes,  (3)  decrease  force  when  both  the  IF  slopes  and  the  IA  slopes  decrease  over  the 
time,  (4)  muscle  fatigue  when  there  is  a  decrease  in  the  IF  slopes  and  an  increase  in  the  I A  slopes. 
A  common  procedure  in  the  joint  analysis  of  frequency  and  amplitude  is  to  count  the  number  of 
points  in  each  stage  (increase  force,  recovery,  decrease  force  and  fatigue)  [20],  [19].  However,  in 
this  research,  instead  of  counting  the  number  of  points,  it  is  proposed  in  each  stage  to  add  the 
square  of  the  distance  of  the  points  to  the  origin  of  the  IA  x  IF  plane  as  a  measure  of  the  energy 
of  the  point  in  that  stage.  The  aim  of  this  procedure  is  to  determine  whether  the  level  of  fatigue 
reported  by  the  subject  is  directly  related  to  the  rates  of  IF  decrease  and  IA  increase. 
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Figure  2.2  Classification  of  the  IF  and  IA  slopes  using  a  joint  analysis  of  frequency  and  amplitude. 
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3 

Results 


This  chapter  describes  the  experimental  results  obtained  by  applying  the  filter  bank  decomposition 
proposed  in  this  research  to  different  SEMG  signals  provided  by  the  AFRL.  These  results  are 
compared  with  those  given  by  the  standard  techniques,  namely  the  spectrogram  and  the  SPWVD. 
Section  3.1  compares  the  IF  and  IA  slopes  estimated  using  the  filter  bank  decomposition  with  the 
slopes  given  by  the  spectrogram  and  the  SPWVD;  this  comparison  is  based  on  the  computation  of 
Pearson’s  correlation  coefficient  and  the  relative  quadratic  difference.  Note  that  the  IF  and  IA 
slopes  are  the  indices  used  to  assess  muscle  fatigue,  which  is  the  object  of  study  in  this  research. 
Using  these  slopes,  the  method  of  joint  analysis  of  frequency  and  amplitude  classifies  the  muscle 
activity  in  each  hour  into  the  following  stages:  increase  force,  recovery,  decrease  force,  or  fatigue. 
The  results  of  these  classifications  are  presented  in  Section  3.2. 

Finally,  the  energy  of  the  IF  and  IA  slopes  for  each  stage  is  proposed  as  a  measure  of  the  intensity 
level  of  increase  force,  recovery,  decrease  force,  or  fatigue.  The  correlation  between  the  measured 
energy  in  the  fatigue  stage  and  the  perceived  levels  of  discomfort  reported  by  the  subjects  is  then 
evaluated  in  Section  3.2. 


3.1  Comparison  of  IF  and  I A  Slopes 

This  section  compares  the  IF  and  IA  slopes  as  computed  by  the  proposed  filter  bank 
decomposition  with  those  given  by  the  spectrogram  and  the  SPWVD  for  three  kinds  of  muscle: 
splenius  capitis,  sterno,  and  trapezius  (see  Fig.  3.4).  As  mentioned  in  Section  2.3,  the  method 
developed  in  this  project  has  the  advantage  of  being  computationally  faster  and  improves  cross 
terms  supression.  However,  it  is  still  important  to  verify  if  the  obtained  slopes  have  the  same  trend 
as  those  given  by  the  standard  spectrogram  and  SPWVD  techniques  as  the  following  examples 
show. 

Example  1  below  illustrates  the  extraction  of  the  IF  and  IA  slopes  of  a  typical  SEMG  signal 
recorded  from  the  right  trapezius  after  hour  2  when  the  subject  was  performing  an  exertion  of 
70%  MVC.  Figure  3.1(a)  shows  the  signal  samples  corresponding  to  3  minutes  at  the  end  of  hour 
2,  while  Fig.  3.1(b)  shows  the  corresponding  IF  and  IA  slopes  obtained  using  the  three  methods 
(spectrogram,  SPWVD,  and  filter  bank  decomposition).  Observe  in  the  top  of  Fig.  3.1(b)  that  the 
IF  slopes  computed  using  the  three  methods  are,  most  of  the  time,  very  close.  Although  the  IA 
slopes,  shown  in  the  bottom  of  Fig.  3.1(b),  are  clearly  different  in  the  peak  region,  the  three 
methods  provide  the  same  trend  of  slopes,  meaning  that  when  the  results  of  one  of  the  methods  is 
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increasing  or  decreasing  with  time  the  same  happens  to  the  other  two  methods.  In  this  way,  in 
terms  of  providing  valid  indices  for  fatigue  assessment,  the  three  methods  can  be  considered 
equivalent  in  this  example,  but  the  proposed  filter  bank  decomposition  method  has  the  advantage 
of  being  faster  and  window-independent. 

A  second  example  is  shown  in  Fig.  3.2.  Figure  3.2(a)  shows  the  signal  samples  extracted  from  the 
right  splenius  during  the  last  3  minutes  of  hour  6  when  the  subject  was  performing  70%  MVC. 
The  corresponding  IF  and  IA  slopes  are  shown  in  Fig.  3.2(b).  As  observed  in  Example  1 
illustrated  in  Fig.  3.1,  the  IF  slopes  obtained  from  the  filter  bank  decomposition  are  equivalent  to 
those  provided  by  the  spectrogram  and  the  SPWVD,  while  the  IA  slopes,  which  differ  in  the  peak 
regions,  present  always  the  same  trends.  Again  it  is  noted  that  the  IF  and  IA  slopes  given  by  the 
proposed  method  are  satisfactory  for  the  assessment  of  muscle  fatigue  with  less  computational 
time  and  with  window-independence. 

The  behavior  verified  in  Figs.  3.1  and  3.2  for  the  proposed  method,  characterized  by  IA  and  IF 
slopes  trends  close  to  those  given  by  the  spectrogram  and  the  SPWVD,  was  typical  of  all  the 
simulations,  provided  that  the  analyzed  signal  does  not  have  high  and  localized  peaks.  These 
peaks  appear  in  some  cases,  probably  due  to  poor  placement  of  the  electrodes  or  bad  contact  of 
the  electrodes  with  the  subject’s  skin.  The  consequences  of  these  wrong  signal  acquisitions  are 
illustrated  in  Example  3,  in  Fig.  3.3.  First,  Fig.  3.3(a)  shows  the  samples  from  the  signal  acquired 
from  the  left  splenius  at  the  3  minutes  at  the  end  of  the  hour  7.  Note  the  high  peaks  around 
minutes  0.3,  1.6,  and  2.4.  In  this  case,  as  it  is  shown  in  the  top  of  Fig.  3.3(b),  the  IF  slopes  differ 
for  all  the  three  evaluated  techniques  (spectrogram,  SPWVD,  and  filter  bank  decomposition).  The 
IA  slopes,  shown  in  the  bottom  of  Fig.  3.3(b),  are  still  very  similar  in  trend  (for  all  three  methods, 
they  increase  or  decrease  in  the  same  time  regions).  However,  the  conclusion  is  that  the  peaks 
observed  in  Fig.  3.3(a)  make  the  IF  slopes  inappropriate  as  indices  of  muscle  fatigue,  as  they  are 
due  to  wrong  measurements. 
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Example  1 
Subject  information 
Subj  ID:  NF  F-l 
Test  #  :  NFF0002 
Cell:  B 

Muscle:  Right  Trapezius 


H2  -  Right  Trapezius 


Figure  3.1  Correct  SEMG  signal  and  correct  estimation  of  its  IF  and  IA.  (a)  SEMG  signal 
recorded  from  the  right  trapezius  at  the  end  of  the  second  hour,  (b)  Computed  slopes  from  the 
IF  and  the  IA  using  three  different  methods  -  spectrogram  (dashed  red  line),  SPWVD  (dotted 
black  line),  and  filter  bank  (continuous  blue  line).  The  IF  and  IA  estimates  from  the  three  methods 
show  similar  trends. 
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Example  2 
Subject  information 
Subj  ID:  NF  F-l 
Test  #  :  NFF0060 
Cell:  E 

Muscle:  Right  Splenius 


H6  -  Right  Splenius 


Time  (min) 

(b) 

Figure  3.2  Correct  SEMG  signal  and  correct  estimation  of  its  IF  and  IA.  (a)  SEMG  signal 
recorded  from  the  right  splenius  at  the  end  of  the  sixth  hour,  (b)  Computed  slopes  from  the  IF 
and  the  IA  using  three  different  methods  -  spectrogram  (dashed  red  line),  SPWVD  (dotted  black 
line),  and  filter  bank  (continuous  blue  line).  The  IF  and  IA  estimates  from  the  three  methods  show 
similar  trends. 
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Example  3 
Subject  information 
Subj  ID:  NF  F-l 
Test  #  :  NFF00060 
Cell:  E 

Muscle:  Left  Splenius 


H7  -  Left  Splenius 


- 

: 

- 

0  0.5  1  1.5  2  2.5  3 

Time  (min) 


(a) 


Figure  3.3  Wrong  SEMG  signal  and  wrong  estimation  of  its  IF  and  IA.  (a)  SEMG  signal  recorded 
from  the  left  splenius  at  the  end  of  the  seventh  hour.  This  signal  shows  high  peaks  and  small 
amplitude  values  at  certain  time.  This  noise  seems  to  be  caused  by  a  bad  contact  of  the  electrodes 
with  the  skin  of  the  subject  or  with  the  EMG  equipment,  (b)  Computed  slopes  from  the  IF  and  the 
IA  using  three  different  methods  -  spectrogram  (dashed  red  line),  SPWVD  (dotted  black  line),  and 
filter  bank  (continuous  blue  line).  The  IF  estimated  from  the  three  methods  show  different  values. 


43 


Two  figures  of  merit  were  used  to  quantify  the  comparison  between  the  slopes  computed  from  the 
subbands  of  the  filter  bank  with  the  slopes  computed  from  either  the  spectrogram  or  SPWVD. 

The  first  figure  of  merit  was  the  relative  quadratic  difference  defined  as 


En=Mb[n\  -  Stf[n ])2 

E»=i(s/ftM)2 


x  100, 


(3.1) 


where  sy&  [n]  is  either  the  IF  slope  or  the  IA  slope  at  time  n  estimated  from  the  filter  bank,  and 
Stf[n]  is  either  the  IF  slope  or  the  IA  slope  at  time  n  estimated  from  either  the  spectrogram  or  the 
SPWVD. 

The  second  figure  of  merit  was  the  Pearson’s  correlation  coefficients  defined  as 


P  =  7 - E  ( sfbH  -  Sfb )  {stf[n\  -  Stf ) , 

L  Vfb  <hf  n=  1 


(3.2) 


where  ay 5,  and  Sfb  are  the  standard  deviations  and  the  sample  means,  respectively,  of  either  the  IF 
or  IA  slopes  computed  from  the  filter  bank.  The  parameters  atf,  and  S/  j  are  the  standard  deviations 
and  the  sample  means,  respectively,  of  either  the  IF  or  IA  slopes  computed  from  either  the 
spectrogram  or  the  SPWVD.  The  relative  quadratic  difference  is  computed  for  n  —  0, 1,  2, . . . ,  L, 
where  L  is  the  number  of  slopes  that  were  computed  over  an  interval  of  the  estimated  IF  or  IA. 
The  results  for  the  relative  quadratic  difference  computed  for  3480  SEMG  signals  provided  by  the 
AFRL  are  summarized  in  Table  3.1.  In  the  case  of  the  IF  slopes  computed  with  the  filter  bank, 
96%  had  differences  not  greater  than  0.5%  with  respect  to  those  slopes  given  by  the  spectrogram, 
and  90%  had  differences  not  greater  than  0.5%  with  respect  to  those  slopes  given  by  the  SPWVD. 
In  the  case  of  the  IA  slopes  computed  with  the  filter  bank,  99%  had  differences  not  greater  than 
0.5%  with  respect  to  those  slopes  given  by  the  spectrogram,  and  100%  had  differences  not  greater 
than  0.5%  with  respect  to  those  slopes  given  by  the  SPWVD. 

The  results  for  the  Pearson’s  correlation  coefficients  computed  for  3480  SEMG  signals  provided 
by  the  AFRL  are  summarized  in  Table  3.2.  Note  that  a  correlation  coefficient  of  1  means  that  the 
slopes  derived  from  the  spectrogram  and  the  SPWVD  follow  the  same  trend  as  the  slopes  derived 
from  the  filter  bank.  In  the  case  of  the  IF  slopes  computed  with  the  filter  bank,  94%  had 
correlation  coefficients  greater  than  0.9  with  respect  to  those  slopes  given  by  the  spectrogram,  and 
88%  had  correlation  coefficients  greater  than  0.9  with  respect  to  those  slopes  given  by  the 
SPWVD.  In  the  case  of  the  IA  slopes  computed  with  the  filter  bank,  97%  had  correlation 
coefficients  greater  than  0.9  with  respect  to  those  slopes  given  by  the  spectrogram,  and  100%  had 
correlation  coefficients  greater  than  0.9  with  respect  to  those  slopes  given  by  the  SPWVD. 

Cohen  has  limited  the  use  of  wavelets  and  filter  bank  in  the  TF  analysis  since  physical  quantities 
derived  from  the  TF  distribution  such  as  the  moments  are  meaningless  or  sometimes  they  do  not 
exist  [7].  However,  from  the  results  of  the  filter  bank  proposed  in  this  research,  it  is  inferred  that 
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slopes  derived  from  the  spectrogram  and  the  SPWVD  show  not  only  identical  trends  but  also 
similar  values  to  those  slopes  derived  from  the  filter  bank.  For  a  few  number  of  slopes,  3%  of  the 
results,  neither  the  filter  bank,  the  spectrogram,  or  the  SPWVD  show  equivalent  slopes.  This  is 
because  a  wrong  acquisition  of  the  signal  causes  high  localized  peaks  in  some  samples,  which  can 
lead  to  invalid  indices  mainly  in  the  case  of  the  IF,  as  shown  in  Fig.  3.3. 


(a)  (b)  (c) 

Figure  3.4  Neck  muscles  involved  in  the  analysis  of  muscle  fatigue,  (a)  Muscle  -  Splenius 
capitis;  Action  -  Extends  and  rotates  cervical  spine,  (b)  Muscle  -  Sterno;  Action  -  Flexes  and 
laterally  rotates  cervical  spine.  Extends  neck  when  neck  is  already  partially  extended;  (c)  Muscle 
-  Trapezius.  Action  -  Stabilizes,  elevates,  retracts,  and  rotates  scapula  (source  -  AFRL). 


Table  3.1  Relative  quadratic  difference  a  of  the  IF  and  IA  slopes  obtained  from  the  spectrogram, 
SPWVD,  and  filter  bank  during  the  analysis  of  3480  SEMG  signals. 


Method 

%  of  IF  slopes 

%  of  IA  slopes 

with  a  <  0.5% 

with  a  <  0.5% 

Spectrogram  x  Filter  Bank 

96 

99 

SPWVD  x  Filter  Bank 

90 

100 

Table  3.2  Pearson’s  correlation  coefficient  p  of  the  IF  and  IA  slopes  obtained  from  the  spectro¬ 
gram,  SPWVD,  and  filter  bank  during  the  analysis  of  3480  SEMG  signals. 


Method 

%  of  IF  slopes 

%  of  IA  slopes 

with  p  >  0.9 

with  p  >  0.9 

Spectrogram  x  Filter  Bank 

94 

97 

SPWVD  x  Filter  Bank 

88 

100 
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3.2  Joint  Analysis  of  Frequency  and  Amplitude 


This  section  describes  the  results  obtained  for  the  classification  of  the  muscle  activity  as  stages  of 
increase  force,  recovery,  decrease  force,  or  fatigue.  This  classification  is  based  on  the  IF  and  IA 
slopes,  which  were  obtained  for  a  subject  during  the  eight-hour  session.  Also,  the  energy  of  the  IF 
and  IA  slopes  is  proposed  as  a  measure  of  the  intensity  level  of  muscle  activity  at  each  stage.  To 
evaluate  the  effectiveness  of  this  measure,  the  correlation  between  the  measured  energy  in  the 
fatigue  stage  and  the  perceived  levels  of  discomfort  reported  by  the  subjects  is  evaluated. 

The  results  of  the  classification  of  the  muscle  activity  stage  for  a  female  subject  performing  an 
exertion  of  70%  MVC  while  using  different  weighted  flight  helmets  A  (3.0  Lb  baseline  CG),  B 
(4.5  Lb  near  head  CG),  C  (4.5  Lb  forward  head  CG),  and  E  (6.0  Lb  forward  head  CG)  are  plotted 
respectively  in  Figs.  3.5,  3.7,  3.9,  and  3.11  (helmet  D  with  4.5  Lb  near  head  CG  was  not 
considered  because  there  was  no  signal  to  analyze).  In  each  figure,  items  (a)  to  (g)  corresponding 
to  each  hour,  present  the  classified  slopes  derived  from  the  filter  bank  decomposition  in  an  IA  x 
IF  plane,  in  which  each  quadrant  corresponds  to  a  different  stage  (increase  force,  recovery, 
decrease  force,  and  fatigue).  In  this  way,  the  muscle  activity  is  classified  as:  (1)  increase  force 
when  both  the  IF  slopes  and  the  IA  slopes  increase  over  time,  (2)  recovery  from  previous  muscle 
fatigue  when  there  is  an  increase  in  the  IF  slopes  and  a  decrease  in  the  IA  slopes,  (3)  decrease 
force  when  both  the  IF  slopes  and  the  IA  slopes  decrease  over  time,  (4)  muscle  fatigue  when  there 
is  a  decrease  in  the  IF  slopes  and  an  increase  in  the  IA  slopes. 

Also  shown  in  Figs.  3.5,  3.7,  3.9,  and  3.11  are  the  energy  levels  (in  blue)  that  were  computed 
from  the  slopes  given  by  the  filter  bank  decomposition.  The  correlation  between  these  energy 
levels  and  the  perceived  levels  of  discomfort  reported  by  the  subjects  is  evaluated  next. 

It  is  inferred  from  Fig.  3.6  that  in  the  fifth  and  sixth  hour  most  of  the  energy  is  concentrated  in 
muscle  fatigue  (28.72%  in  the  fifth  hour  and  25.17%  in  the  sixth  hour);  however,  the  first  and 
second  hour  show  less  energy  in  the  stage  fatigue  (3.55%  in  the  first  hour  and  13%  in  the  second 
hour).  According  to  these  results,  the  subject  is  expected  to  report  more  fatigue  in  the  sixth  hour 
than  in  the  second  hour.  The  previous  assumption  is  confirmed  with  the  perceived  levels  of 
discomfort  reported  by  the  subject  when  wearing  helmet  A.  In  the  questionnaire  the  subject 
reported  for  the  following  regions  (head,  upper  neck,  lower  neck,  shoulders,  upper  back,  lower 
back)  in  the  second  hour:  (1)  no  discomfort  in  the  upper  back,  upper  neck,  and  lower  neck;  (2)  no 
hot  spots  and  no  numbness  in  any  region.  In  addition,  the  subject  reported  no  headache,  relaxed 
mental  state,  alert  mental  frame  of  mind,  and  attentive  concentration  level.  In  contrast,  during  the 
sixth  hour  the  subject  reported  soreness  and  moderate  hot  spots  in  all  the  regions,  but  no 
numbness  in  any  region.  In  addition,  the  subject  reported  minimal  headache,  slightly  tense,  tired 
mental  frame  of  mind,  and  distracted  concentration  level. 

Figure  3.8  shows  the  analysis  of  muscle  fatigue  when  the  subject  wore  helmet  B.  In  this  analysis, 
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the  right  splenius  and  right  trapezius  were  the  muscles  most  fatigued  in  the  first,  second,  and 
fourth  hour.  In  the  sixth  hour  the  right  sterno  is  the  only  muscle  showing  a  recovery  from  previous 
muscle  fatigue.  Most  of  the  neck  muscles  present  fatigue  at  the  end  of  the  fifth,  sixth,  and  seventh 
hour.  In  the  questionnaire  the  subject  reported  for  the  folowing  regions  (head,  upper  neck,  lower 
neck,  shoulders,  upper  back,  lower  back)  in  the  second  and  fourth  hour:  (1)  no  discomfort  in  the 
head,  upper  neck,  lower  neck,  shoulders,  upper  back,  and  lower  back;  (2)  no  hot  spots  and  no 
numbness  in  any  region.  In  addition,  the  subject  reported  no  headache,  relaxed  mental  state,  alert 
mental  frame  of  mind,  and  attentive  concentration  level.  In  contrast,  during  the  sixth  hour  the 
subject  reported:  (1)  soreness  in  the  head,  upper  neck,  lower  neck,  shoulders  and  upper  back,  and 
discomfort  in  the  lower  back;  (2)  moderate  hot  spots  in  all  the  regions,  but  no  numbness  in  any 
region.  In  addition,  the  subject  reported  severe  headache,  slightly  tense,  tired  mental  frame  of 
mind,  and  distracted  concentration  level. 

In  the  case  of  helmet  C  ,  the  results  derived  from  the  filter  bank  could  not  be  correlated  to  the 
perceived  levels  of  discomfort  reported  by  the  subject  because  these  discomfort  levels  were  not 
found  in  the  data  provided  by  the  AFRL.  However,  from  the  analysis  of  helmet  C,  it  can  be 
inferred  that  the  right  splenius  is,  during  most  of  the  time  except  for  the  fourth  hour,  the  most 
fatigued  muscle.  In  addition,  it  can  be  concluded  that  the  neck  muscles  in  the  seventh  hour  are 
more  fatigued  than  in  the  previous  hours. 

For  the  analysis  of  muscle  fatigue  when  the  subject  wore  helmet  E,  it  can  be  inferred  from  Fig. 
3.12  that  most  of  the  energy  related  to  fatigue  is  concentrated  in  the  fifth  hour  and  sixth  hour 
(32.8%  in  the  fifth  hour  and  32.5%  in  the  sixth  hour).  The  second  and  the  fourth  hour  show  less 
fatigue  (8%  in  the  second  hour  and  1 1%  in  the  fourth  hour).  The  right  splenius  is  the  most 
fatigued  muscle  during  the  fifth  and  sixth  hour.  In  regard  to  the  questionnaire,  the  subject  reported 
for  the  following  regions  (head,  upper  neck,  lower  neck,  shoulders,  upper  back,  lower  back)  in  the 
second  hour:  (1)  discomfort  in  the  head,  upper  neck,  lower  neck,  shoulders  and  upper  back,  and 
no  discomfort  in  the  lower  back;  (2)  moderate  hot  spots  in  the  upper  neck,  lower  neck,  shoulders 
and  upper  back;  (3)  no  numbness  in  any  region.  In  addition,  the  subject  reported  minimal 
headache,  slightly  tense,  tired  mental  frame  of  mind,  and  attentive  concentration  level.  In  the 
fourth  hour  the  subject  reported:  (1)  soreness  in  head,  upper  neck,  lower  neck,  shoulders  and 
upper  back,  and  no  discomfort  in  the  lower  back;  (2)  moderate  hot  spots  in  the  head,  upper  neck, 
lower  neck,  shoulders  and  upper  back,  and  no  hot  spots  in  the  lower  back;  (3)  no  numbness  in  any 
region.  In  addition,  the  subject  reported  severe  headache,  slightly  tense,  tired  mental  frame  of 
mind,  and  distracted  concentration  level.  Finally,  in  the  sixth  hour  the  subject  felt  a  stronger 
fatigue  reporting:  (1)  pain  in  the  head,  upper  neck,  lower  neck,  shoulders  and  upper  back,  and 
discomfort  in  the  lower  back;  (2)  severe  hot  spots  in  the  head,  upper  neck,  lower  neck,  shoulders 
and  upper  back,  and  moderate  hot  spots  in  the  lower  back;  (3)  no  numbness  in  any  region.  In 
addition,  the  subject  reported  severe  headache,  restless  mental  state,  exhausted  mental  frame  of 
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mind,  and  loss  of  focus  in  the  concentration  level.  Pain  reported  by  this  subject  correlates  the 
results  obtained  from  the  filter  bank. 

During  this  analysis  of  muscle  fatigue,  it  is  noted  that  the  energy  levels  of  fatigue  are  not 
consistent  for  all  the  helmets.  It  was  expected  that  a  higher  energy  in  the  stage  of  fatigue  would  be 
obtained  when  the  subject  wore  helmet  E  instead  of  A.  In  addition,  if  a  set  of  muscles  are  fatigued 
in  the  fifth  and  the  sixth  hour,  it  was  expected  that  the  same  set  of  muscles  would  be  fatigued  in 
hour  7.  However,  this  did  not  occur  in  the  analysis. 
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Figure  3.5  Joint  analysis  of  frequency  and  am¬ 
plitude  for  subject  ID  #  NF  F-l,  test  #  NFF0001, 
cell  A  (3.0  Lb  baseline  CG).  Stages  (I)  muscle  in¬ 
crease  force,  (II)  recovery,  (III)  muscle  decrease 
force,  (IV)  fatigue,  (a)  hour  1,  (b)  hour  2,  (c)  hour 
3,  (d)  hour  4,  (e)  hour  5,  (f)  hour  6,  (g)  hour  7. 
The  numbers  inside  each  quadrant  (in  blue)  cor¬ 
respond  to  the  energy  of  the  slopes  in  that  stage. 
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eight  hour  session. 
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Figure  3.7  Joint  analysis  of  frequency  and  am¬ 
plitude  for  subject  ID  #  NF  F-l,  test  #  NFF0002, 
cell  B  (4.5  Lb  near  head  CG).  Stages  (I)  mus¬ 
cle  increase  force,  (II)  recovery,  (III)  muscle  de¬ 
crease  force,  (IV)  fatigue,  (a)  hour  1,  (b)  hour 
2,  (c)  hour  3,  (d)  hour  4,  (e)  hour  5,  (f)  hour  6, 
(g)  hour  7.  The  numbers  inside  each  quadrant  (in 
blue)  correspond  to  the  energy  of  the  slopes  in 
that  stage. 
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previous  muscle  fatigue  during  the  sixth  hour.  Most  of  the  neck  muscles  present  fatigue  at  the  end  of  the  fifth,  sixth  and  seventh  hour. 
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Figure  3.9  Joint  analysis  of  frequency  and  am¬ 
plitude  for  subject  ID  #  NF  F-l,  test  #  NFF0024, 
cell  C  (4.5  Lb  forward  head  CG).  Stages  (I)  mus¬ 
cle  increase  force,  (II)  recovery,  (III)  muscle  de¬ 
crease  force,  (IV)  fatigue,  (a)  hour  1,  (b)  hour 
2,  (c)  hour  3,  (d)  hour  4,  (e)  hour  5,  (f)  hour  6, 
(g)  hour  7.  The  numbers  inside  each  quadrant  (in 
blue)  correspond  to  the  energy  of  the  slopes  in 
that  stage. 
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Figure  3.11  Joint  analysis  of  frequency  and  am¬ 
plitude  for  subject  ID  #  NF  F-l,  test  #  NFF0060, 
cell  E  (6.0  Lb  forward  head  CG).  Stages  (I)  mus¬ 
cle  increase  force,  (II)  recovery,  (III)  muscle  de¬ 
crease  force,  (IV)  fatigue,  (a)  hour  1,  (b)  hour 
2,  (c)  hour  3,  (d)  hour  4,  (e)  hour  5,  (f)  hour  6, 
(g)  hour  7.  The  numbers  inside  each  quadrant  (in 
blue)  correspond  to  the  energy  of  the  slopes  in 
that  stage. 
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hour. 


4 

Conclusion 


This  project  proposes  a  technique  for  estimating  the  instantaneous  frequency  (IF)  and  the 
instantaneous  amplitude  (IA)  of  a  surface  electromyographic  (SEMG)  signal  by  decomposing  it 
using  a  32-channel  filter  bank.  Both  the  IF  and  IA  should  be  used  for  the  assessment  of  muscle 
fatigue,  since  in  muscle  fatigue  the  SEMG  signal  suffers  a  decrease  in  the  IF  and  an  increase  in 
the  IA.  Then,  by  computing  the  IF  and  IA  slopes,  which  measure  their  decrease  or  increase  rate 
over  time,  the  method  of  joint  analysis  of  frequency  and  amplitude  is  used  to  classify  the  muscle 
activity  into  one  of  the  following  stages:  increase  force,  recovery,  decrease  force,  or  fatigue. 

The  proposed  method  for  the  assessment  of  muscle  fatigue  has  the  advantage  of  being 
computationally  less  expensive  than  other  standard  techniques,  such  as  the  Wigner-Ville 
distribution  (WVD),  the  spectrogram,  the  Choi-Williams  distribution  (CWD),  and  the  smoothed 
pseudo  Wigner-Ville  distribution  (SPWVD).  The  filter  bank  is  window  independent  and  provides 
a  better  cross  term  supression  than  the  other  evaluated  techniques.  In  addition,  the  filter  bank  can 
be  implemented  by  iterating  a  two  channel  filter  bank  used  to  compute  the  wavelet  coefficients  of 
the  SEMG  signal  in  a  scheme  also  known  as  wavelet  packets.  The  advantage  of  using  a  wavelet 
decomposition  is  that  the  estimation  of  the  IF  and  IA  can  be  combined  with  other  techniques 
based  on  wavelet  transforms  such  as  baseline  drift  removal  and  signal  noise  reduction. 

In  the  proposed  method,  the  IF  is  estimated  from  the  subbands  of  the  filter  bank  as  the  first  order 
moment  of  the  time  subband  representation  at  a  certain  point  in  time.  The  IA  is  estimated  also 
from  the  subbands  of  the  filter  bank,  according  to  a  proposed  technique  based  on  applying  the 
outputs  of  the  filters  to  successive,  overlapping  adaptive  windows  and  determining  the  maximum 
coefficient  from  each  window.  The  sum  of  the  maximum  coefficients  for  the  different  subbands 
then  corresponds  to  a  single  IA  estimate.  From  the  computed  IF  and  IA  estimates,  the  slopes 
needed  for  the  classification  of  muscle  activity  are  then  obtained  by  using  a  first  order  linear 
regression  model  for  intervals  of  one  minute  every  one  second. 

For  evaluating  the  IF  and  IA  slopes  computed  from  the  filter  bank  decomposition,  two  other 
standard  techniques,  the  spectrogram  and  the  SPWVD,  were  also  implemented  since  they  were 
shown  to  provide  a  better  estimation  of  the  IF  and  IA  than  the  WVD  and  the  CWD.  From  the 
spectrogram  and  the  SPWVD,  the  IF  was  also  estimated  as  the  first  order  moment  of  the  TF 
distribution,  and  the  IA  was  estimated  as  the  squared  root  of  the  time  marginal. 

During  the  performance  tests,  the  IF  and  IA  slopes  were  computed  from  the  filter  bank 
decomposition,  the  spectrogram,  and  the  SPWVD  of  3480  SEMG  signals  provided  by  the  AFRF1 . 

'The  Kynesiology  Lab  from  UTEP  recorded  SEMG  signals  from  six  lower  back  muscles,  4  hip  muscles,  and  4 
hamstring  muscles.  The  physical  routines  performed  by  the  subjects  in  the  Kynesiology  Lab  to  analyze  muscle  fatigue 
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Slopes  obtained  from  the  filter  bank  decomposition  were  correlated  and  compared  to  the  slopes 
obtained  from  the  spectrogram  and  the  SPWVD.  In  the  case  of  the  IF  slopes  computed  with  the 
filter  bank,  94%  showed  correlation  coefficients  greater  than  0.9  with  respect  to  those  slopes  given 
by  the  spectrogram,  and  88%  showed  correlation  coefficients  greater  than  0.9  with  respect  to 
those  slopes  given  by  the  SPWVD.  In  the  case  of  the  IA  slopes  computed  with  the  filter  bank, 

97%  showed  correlation  coefficients  greater  than  0.9  with  respect  to  those  slopes  given  by  the 
spectrogram,  and  100%  showed  correlation  coefficients  greater  than  0.9  with  respect  to  those 
slopes  given  by  the  SPWVD.  These  results  verify  that  the  slopes  obtained  from  the  filter  bank 
show  a  strong  correlation  and  the  same  trend  as  those  given  by  the  standard  spectrogram  and  the 
SPWVD  techniques. 

The  relative  quadratic  differences,  when  comparing  the  IF  and  IA  slopes  computed  from  the  filter 
bank  with  those  slopes  computed  from  the  spectrogram  and  the  SPWVD,  showed  that  89.6%  of 
the  differences  were  less  than  0.5%.  This  means  that  the  slopes  obtained  from  the  filter  bank  not 
only  show  the  same  trend  but  also  close  values  to  those  given  by  the  spectrogram  and  the 
SPWVD. 

For  a  few  number  of  slopes,  3%  of  the  results  obtained  from  the  data  provided  by  the  AFRL, 
neither  the  filter  bank,  the  spectrogram,  or  the  SPWVD  showed  equivalent  slopes.  This  was 
attributed  to  the  fact  that  the  corresponding  SEMG  signals  show  high  peak  values  at  some 
samples  due  to  poor  placement  of  the  electrodes  or  bad  contact  of  the  electrodes  with  the  skin  of 
the  subject. 

By  using  a  joint  analysis  of  frequency  and  amplitude,  the  IF  and  IA  slopes  obtained  from  the  filter 
bank  decomposition  are  classified  into  one  of  the  following  four  stages:  increase  force,  recovery, 
muscle  decrease  force  or  fatigue.  The  energy  of  the  IF  and  IA  slopes  for  each  stage  is  proposed  as 
an  index  to  measure  the  intensity  level  of  increase  force,  recovery,  decrease  force,  or  fatigue.  The 
correlation  between  the  measured  energy  in  the  fatigue  stage  and  the  perceived  levels  of 
discomfort  reported  by  the  subjects  were  then  evaluated. 

An  extensive  evaluation  of  the  3480  SEMG  signals  showed  that  the  proposed  filter  bank 
decomposition  is  a  suitable  tool  for  determining  muscle  fatigue.  In  all  the  cases  that  the  subject 
reported  any  discomfort  or  fatigue  for  a  group  of  muscles,  the  proposed  filter  bank  provided  an 
index  indicating  muscle  fatigue.  However,  this  index  is  not  linearly  proportional  to  the  intensity 
of  fatigue  felt  by  the  subject.  For  instance,  muscle  fatigue  was  analyzed  for  a  subject  performing 
70%  MVC  and  wearing  helmets  A,  B,  C,  and  E  during  an  eight-hour  session.  In  this  analysis  was 
shown  for  all  the  analyzed  helmets  (A,  B,  C,  and  E)  that  the  proposed  filter  bank  decomposition 
always  identified  muscle  fatigue  which  matched  the  reported  perceived  levels  of  discomfort. 

were  much  lighter  than  the  flight  routines  performed  by  the  subjects  in  the  AFRL  experiment.  As  a  result,  perceived 
levels  of  discomfort  reported  by  the  subjects  in  the  UTEP’s  Kinesiology  Lab  experiment  reported  no  cases  of  muscle 
fatigue,  in  contrast  to  the  AFRL  experiment  that  reported  cases  of  muscle  fatigue.  Therefore,  the  assessment  of  muscle 
fatigue  used  in  this  research  was  applied  only  to  SEMG  signals  provided  by  the  AFRL. 
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Higher  indices  of  muscle  fatigue  were  expected  when  the  subject  wore  helmet  E  compared  to 
helmet  A,  since  helmet  E  was  heavier.  However,  the  indices  of  muscle  fatigue  obtained  for  helmet 
E  were  lower  than  helmet  A. 

A  limitation  of  the  experimental  results  is  that  the  levels  of  discomfort  were  reported  for  certain 
regions  (head,  upper  neck,  lower  neck,  shoulders,  upper  back,  and  lower  back)  related  to  the 
analyzed  neck  muscles  but  not  specifically  for  each  muscle.  Moreover,  the  discomfort  levels  are 
characteristic  of  the  subject  and  can  change  depending  on  the  mood,  time,  weather,  and  other 
factors  that  are  involved  during  the  collection  of  the  SEMG  signals.  Therefore,  it  is  important  to 
correlate  SEMG  signals  to  other  biomedical  measurements  such  as  blood  sample,  heart  rate,  and 
electroencephalographic  (EEG)  signals  in  order  to  create  a  more  reliable  index  to  measure  the 
intensity  of  muscle  fatigue  felt  by  the  subject. 
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Appendix  A 
Toolbox 
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Algorithm  1  Compute  the  instantaneous  frequency  and  the  instantaneous  amplitude 


Input: 

-  Define  the  signal  x[n],\/  n  =  0, 1, 2,  3, . . . ,  L  —  1,  where  L  is  the  length  of  the  signal 

-  Define  a  filter  bank  h[f,  n],  V  n  —  0, 1, 2,  3, . . . ,  M  —  1  and  /  =  0, 1,  2, 3  . . . ,  F  —  1,  where  M 
is  the  length  of  each  filter  and  F  is  the  number  of  subbands.  The  filter  /0  is  a  low  pass  filter 
and  the  filter  /f-i  is  a  high  pass  filter.  In  this  tool  box  all  the  filters  have  the  same  length 

Output: 

-  Instantaneous  frequency  f[n] 

-  Instantaneous  amplitude  ia[n] 

1:  Compute  the  subband  coefficients  y[f,  n]  convolving  the  input  x[n]  with  the  filter  bank  h[f,  n] ; 
y[f,  n]  =  x[n]  *  h[f,  n] .  The  resulting  output  y[f,  n]  is  defined  V  n  —  0, 1,  2,  3, . . . ,  L  +  M  —  2 
and  /  =  [0, 1,  2,  3, . . . ,  F  —  l]^p 
2:  Define  the  output  y[f,  n\  as 

y[f,  n]  =  y[f ,  n  +  f],  V  n  =  0, 1, 2, 3, . . . ,  L  -  1 
3:  Compute  the  time  subband  representation  as 
y[f,n]  =  (: y[f,n ])2 

4:  Compute  the  time  marginal  m[n\  as 

F—l 

m[n]  =  V[f,n\ 

f=o 

5:  Compute  the  instantaneous  frequency  f\n]  as 
?r  n  E/=o  fy[f,n] 
fln]  =  m\a] 

6:  Compute  the  instantaneous  amplitude  ia[n]  as 
7:  for  /  =  0  to  F  —  1  do 
8:  r  =  ceil(-^j-),  k  =  0,  flag=  1 

9:  while  flag^  0  do 

10:  sf)k[n\  =  y[f,  n  +  k}-  V  n  =  0, 1,  2,  3, . . . ,  r  -  1 

11:  V’t/jfc]  =  max(s/ifc[n]) 

12:  r  —  L  +  M  —  1  —  k 

13:  if  r  =  0  then 

14:  flag  =  0 

15:  else 

16:  k  =  k  +  1 

17:  end  if 

18:  end  while 

19:  end  for 

20:  ia[n  =  k]  =  ^  f>[f,  k] 
f= o 
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Matlab  Functions  Content 


The  following  functions  shown  in  Table  A.  1  are  implemented  in  Matlab.  These  functions  are  used 
in  this  research  to  assess  muscle  fatigue. 


Table  A.l  Matlab  functions  used  in  the  assessment  of  muscle  fatigue. 


Matlab  Function 

Purpose 

gfb 

Generate  the  analysis  filter  bank 

f_filter_bank_ana 

Compute  the  subband  outputs  of  the  filter  bank 

tsub 

Arrange  all  the  computed  subband  coefficients  in  a  matrix  form 

inst  freq 

Compute  the  instantaneous  frequency  of  a  given  signal 

inst_amp 

Compute  the  instantaneous  amplitude  of  a  given  signal 

windowing 

Compute  the  start  and  end  indexes  of  the  location  of  a  rectangular  moving  window 

linear  j-egression 

Estimate  the  slopes  and  points  of  interception 

jasa_analysis 

Compute  the  JASA  method  to  an  electromyography  signal 
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gfb 


Purpose 

Generate  the  analysis  filter  bank. 

Synopsis 

[fltana,  indf ltana] =gfb (n_sb) 

Description 

gfb  generates  the  analysis  filterbank.  For  this  toolbox  the  number  of  sub¬ 
bands  n  sb  can  be  16,  32  or  64.  We  can  represent  the  filterbank  in  matrix 
form  by  h[f,  n]=f  ltana(/My  +  71  +  1)  for  /=0, 1,  2,  3, . , . ,  F  —  1  and 
n=0, 1,  2.  3, ....  Mj  —  1  where  F  is  the  number  of  subbands  n_sb  and  Mj  is  the 
length  of  the  filter  in  subband  /.  The  length  of  the  filter  /  can  be  computed  from 
indf  ltana  using  M/=indf  ltana(/  +  1)  — indf  ltana(/)  for  /  >  0,  and 
indf  ltana  (/  +  1)  for  /= 0.  The  output  indf  ltana  is  required  only  when 
the  filters  have  different  length.  In  this  project  all  the  filters  used  had  length 
equal  to  128  coefficients. 


Name  Description  Default  value 

n_sb  number  of  subbands  32 

fltana  column  vector  with  the  coefficients  of  the  analysis 
filter  bank 

indfltana  column  vector  with  the  index  of  the  end  of  each 
filter 
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Example 

[fltana,  indfltana] =gfb (32) ; 
hn=reshape (fltana, 128, 32) ; 

Hw=fft (hn, 512) ; 

w=l inspace  (0,  1,256) ; 

plot (w,abs (Hw  (256:511,  :) ) ) ; 

x label ( ' $\frac( \ omega } { \pi } $ ' ,  'font size', 14, . . . 
' interpreter' ,  ' latex' ) ; 

y label ( ' $ | H ( \ omega) 1$',  'font size',  14,... 

' interpreter' ,  ' latex' ) ;  grid  on 


7 r 

Figure  A.l  Frequency  response  of  the  analysis  filterbank  (32  subbands)  used  in  the  estimation  of 
the  IF  and  the  IA  of  the  SEMG  signal. 
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f_filter_bank_ana 


Purpose 

Compute  the  subband  outputs  of  the  filter  bank. 

Synopsis 

[ sO ,  isO ] =f_f ilter_bank_ana (sig, fit ana, indf It ana) 

Description 

f_f ilter Jbank_ana  computes  the  subband  coefficients  sO  by  convolving 
the  signal  sig  with  the  filter  coefficients  fit  an  a.  Representing  the  input 
signal  as  x[n],  for  n=0, 1,  2, 3, . . , ,  L  —  1  where  L  is  the  length  of  the  signal,  then 
we  have  x[n)=siq(n  +  1).  The  output  y[f,  n],  for  /= 0, 1,  2,  3, . . . ,  F  —  1  and 
n= 0, 1,2,3 +  Mf  —  2  where  F  is  the  number  of  subbands  n_sb  and  M f 
is  the  length  of  the  filter  in  subband  /,  is  given  by  y[f,  n]=s  0  (i  s  0(f)  +  n  +  1), 
for  /  >  0  and  s0(n  +  1),  for  /= 0. 

Name  Description  Default  value 

sig  column  vector  of  size  Lx  1  with  the  coefficients 

of  the  input  signal 

fit  an  a  column  vector  of  size  F  Mf  x  1  with  the  coeffi¬ 

cients  of  the  analysis  filter  bank 

indfltana  column  vector  of  size  Fxl  with  the  index  of  the 
end  of  each  filter 

sO  column  vector  of  size  Y^f=o(L  +  Mf  —  1)  x  1  with 

the  coefficients  of  each  subband 

i  s  0  column  vector  of  size  Fxl  with  the  index  of  the 

end  of  each  subband 


Example 

load  ' NFFO 0 0 l_Hl_7 0_EMG' 

[sig] =data ( : , 2 ) ; 

[fltana,  indfltana] =gfb (32) ; 

[ sO ,  isO ] =f_f ilter_bank_ana (sig, fltana, indfltana) ; 

See  Also 

gfb . 
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tsub 


Purpose 

Arrange  all  the  computed  subband  coefficients  in  a  matrix  form  of  size  F  x  L. 

Synopsis 

[tfr] =tsub  (sO,  isO,  shift); 

Description 

tsub  arranges  in  a  time  subband  matrix  representation  all  the  subband 
coefficients.  This  function  assumes  that  all  the  filters  have  equal  length. 
The  time  subband  coefficients  y[f,  n]  for  /=0, 1,  2,  3, . . . ,  F  —  1  and 
n=0,l,2,3,...,L  +  Mf  —  2  where  F  is  the  number  of  subbands  n_sb, 
L  is  the  length  of  the  signal  sig  and  Mf  is  the  length  of  the  filter  in  subband  / 
is  given  by  y[f,  n]=tf  r(/  +  1,  n  +  1). 

Name  Description  Default  value 

FA 

s  0  column  vector  of  size  (L  +  Mf  —  1)  x  1  with 

/= o 

the  coefficients  of  each  subband 

i  s  0  column  vector  of  size  F  x  1  with  the  index  of  the 

end  of  each  subband 

shift  vector  of  length  two  with  the  number  of  coeffi¬ 
cients  to  remove  from  the  output  s  0  to  compensate 
the  increase  in  the  length  of  the  subbands  produced 
by  the  convolution  operation 
tfr  time  subband  matrix  representation 


Example 

load  ' NFF0001_H1_70_EMG' 

[sig] =data ( : , 2 ) ; 

[fltana,  indfltana] =gfb (32) ; 

N=length (fltana) /length (indfltana) ; 
shift= [ round (N/2 ) , round (N/2 ) ] ; 

[ sO ,  isO ] =f_f ilter_bank_ana (sig, fltana, indfltana) ; 
[tfr] =tsub  (sO,  isO,  shift); 

See  Also 

gfb,  f_f i lter_bank_ana . 
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instJreq 


Purpose 

Compute  the  instantaneous  frequency  of  a  given  signal. 

Synopsis 

[ fm] =inst_f req (t f r ) ; 


Description 

inst  freq  computes  the  instantaneous  frequency  fm  from  the  time  sub¬ 
band  matrix  representation  tfr.  The  instantaneous  frequency  f[n]  at  time 
n=0, 1,  2,  3, . . . ,  L  —  1  is  given  by  f[n]=fm(n  +  1). 


Name 

tfr 

fm 


=  Sfo  fvlf,  n 
J  E/i'sl f,n 


Description 

time  subband  matrix  representation  of  size  F  x  L 
column  vector  of  size  Lx  1  with  the  instantaneous 
frequency  components 


Default  value 
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Example 

load  ' NFFO 0 0 l_Hl_7 0_EMG' 

[ sig] =data ( : , 2 ) ; 

[fltana,  indfltana] =gfb (32) ; 

N=length (fltana) /length (indfltana) ; 
shift= [ round (N/2 ) , round (N/2 ) ] ; 

[ sO ,  isO ] =f_f ilter_bank_ana (sig, fltana, indfltana) ; 
[tfr] =tsub  (sO,  isO,  shift); 

[ fm] =inst_f req  (tf r)  ; 

t= ( 1 : length  (sig) )/(  60*1024 ) ;  plot(t,  fm.*1024); 
axis([l  3  0  511]);  grid  on; 
xlabel('Time  (min)'',  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

ylabel (' Frequency  (Hz)',  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

See  Also 

gfb,  f_f i lter_bank_ana,  tsub. 
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Figure  A.2  Computed  instantaneous  frequency  of  the  input  signal  s  ig  by  using  32  channels  filter 
bank. 
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inst_amp 


Purpose 

Compute  the  instantaneous  amplitude  of  a  given  signal. 

Synopsis 

[ i_amp] =inst_amp (t f r ) ; 


Description 

inst  amp  computes  the  instantaneous  amplitude  i  amp  from  the  time  sub¬ 
band  matrix  representation  tfr.  The  instantaneous  amplitude  ia[n]  at  time 
n=0, 1,  2,  3, . . . ,  L  —  1  is  given  by  ia[n]=i_amp(n  +  1). 


Name 

Description  Default  value 

tfr 

time  subband  matrix  representation  of  size  F  x  L 

i_amp 

column  vector  of  size  Lx  1  with  the  instantaneous 
amplitude  components 
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Example 

load  ' NFFO 0 0 l_Hl_7 0_EMG' 

[ sig] =data ( : , 2 ) ; 

[fltana,  indfltana] =gfb (32) ; 

N=length (fltana) /length (indfltana) ; 
shift= [ round (N/2 ) , round (N/2 ) ] ; 

[ sO ,  isO ] =f_f ilter_bank_ana (sig, fltana, indfltana) ; 
[tfr] =tsub  (sO,  isO,  shift); 

[ i_amp] =inst_amp (t f r ) ; 

t= ( 1 : length (sig) )/( 60*1024 ) ;  plot  (t,  i_amp) ; 
axis([l  3  0  max (i_amp) ] ) ;  grid  on; 
xlabel('Time  (min)',  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

ylabel (' Amplitude' ,  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

See  Also 

gfb,  f_f i lter_bank_ana,  tsub. 
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Figure  A.3  Computed  instantaneous  amplitude  of  the  input  signal  s  ig  by  using  32  channels  filter 
bank. 
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windowing 


Purpose 

Compute  the  start  and  end  indexes  of  the  location  of  a  rectangular  moving  win¬ 
dow  for  a  specified  shift  amount. 


Synopsis 

[ index_start , 
[ index_start , 
[ index_start , 


index_end] =windowing (nc,  L, 
index_end] =windowing (nc,  L, 
index_end] =windowing (nc,  L) 


shift , 
shift ) 


option) 


Description 

windowing  computes  the  start  and  end  indexes  of  the  location  of  a  rectangular 
moving  window  for  a  specified  shift  amount.  This  function  is  useful  when  we 
want  to  apply  a  rectangular  window  to  a  signal  of  length  L. 


Name 

Description 

Default  value 

nc 

L 

option 

size  of  the  rectangular  window,  for  1  <  nc  <  L 
length  of  the  signal 

if  — >  1  (index_end  —  index_start)  is  always 

1 

shift 

nc.  The  rectangular  window  is  not  zero  padded  at 
the  end  of  the  signal. 

if  — >  0  (index.end  —  index.start)  is  not  al¬ 
ways  nc.  The  rectangular  window  is  zero  padded 
at  the  end  of  the  signal, 
shifted  window,  for  1  <  shift  <  nc 

nc 

index_st  art 

index_end 

column  vector  with  the  start  indexes 
column  vector  with  the  end  indexes 
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Example 


%  Let  us  suppose  that  we  want  to  apply  a  rectangular 
%  window  of  5  points  with  a  shift  of  3  to  a  signal  of 
%  length  15 
nc=5;  L=15;  shift=3; 

%  when  option=l 
opt ion=l ; 

[ index_start_l , index_end_l ] =windowing (nc,  L,  shift,  option) 
%  when  option=0 
opt ion=0 ; 

[ index_start_0 , index_end_0 ] =windowing (nc,  L,  shift,  option) 

%  Output  of  the  code  above  when  option=l 
display (' index_start_l ' ) ; 

index_start_l  = 


1  4  7  10 

display ( ' index_end_l ' ) ; 
index_end_l  = 


5  8  11  14 

%  Output  of  the  code  above  when  option=0 
display (' index_start_0 ' ) ; 

index_st art_0  = 

1  4  7  10  13 

display ( ' index_end_0 ' ) ; 

index_end_0  = 

5  8  11  14  15 
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linear  .regression 


Purpose 

Estimate  the  slopes  and  points  of  interception. 


Synopsis 

[p] =linear_regression (s,  Fs,  ns,  shift) 
[p] =linear_regression (s,  Fs,  ns) 


Description 

linear.regression  estimates  the  slopes  (a*)  and  points  of  interception  (&*) 
of  a  piecewise  approximation  att  +  b,,  for  the  time  interval  t  e  [ti0,£ii]  of  the 
input  data. 


Name  Description  Default  value 


s 

Fs 

ns 

shift 

input  data 

sampling  frequency 

size  of  the  interval  [ti0,  tn] 

shifted  window,  for  1  <  shift  <  length  (s)  . 

The  intervals  of  the  linear  approximation  may 
overlap 

1 

P 

matrix  where  the  first  column  corresponds  to  the 
slopes  and  the  second  column  to  the  points  of  in¬ 
terceptions 
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Example 

load  lt_splenius 

%  to  compute  the  instantaneous  frequency  slopes 
s=fm;  %  the  instantaneous  frequency  fm  is  the  output 
%  of  the  function  inst_freq.m 
Fs  =  l 02  4 ; 

ns=60*1024;  %1  sec  =  1024  samples;  then  60*1024  =  1  minute 
shift=1024;  %shift  every  minute 
[p] =linear_regression (s,  Fs,  ns,  shift); 
t=l inspace (l,3,size(p,l) ) ; 

plot  (t,  p  (  : , 1) , '  r-' , ' linewidth' , 3) ;  grid  on 
xlabel('Time  (min)',  '  f ontsize' ,  14,  ... 

' interpreter' ,  ' latex' ) ; 

ylabel  (' Slope  (Hz/min) ' ,  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

%  to  compute  the  instantaneous  amplitude 

s=i_amp;  %  the  instantaneous  amplitude  i_amp  is 

%  the  output  of  the  function  inst_amp.m 

Fs=1024 ; 

ns=60*1024;  %1  sec  =  1024  samples;  then  60*1024  =  1  minute 
shift=1024;  %shift  every  minute 
[p] =linear_regression (s,  Fs,  ns,  shift); 
t=l in space (l,3,size(p,l)); 

plot  (t,  p (:,  1 )* 60*100 ,' r-' ,' linewidth' , 3) ;  grid  on 
xlabel('Time  (min)',  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 

ylabel (' Slope ( \%/min) ' ,  'fontsize',  14,  ... 

' interpreter' ,  ' latex' ) ; 
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Figure  A.4  Instantaneous  frequency  slopes  from  the  estimated  IF. 


Figure  A.5  Instantaneous  amplitude  slopes  from  the  estimated  IA. 
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jasa_analysis 


Purpose 

Compute  a  joint  analysis  of  the  spectrum  and  the  amplitude  (JASA)  of  a  elec¬ 
tromyography  signal. 

Synopsis 

[i_f,  r,  d_f,  f ] = jasa_analysis (s_inst_amp,  s_inst_freq) 

Description 

jasa.analysis  computes  a  joint  analysis  of  the  spectrum  and  the  amplitude 
of  a  electromyography  signal  (JASA)  from  the  instantaneous  frequency  and  the 
instantaneous  amplitude  slopes  to  compute  indexes  of  muscle  increase  force, 
recovery,  muscle  decrease  force  and  fatigue. 


Name  Description  Default  value 

s_inst_amp  instantaneous  amplitude  slopes 
s_inst_f  req  instantaneous  frequency  slopes 
i_f  Index  corresponding  to  the  energy  of  muscle  in¬ 

crease  force 

r  Index  corresponding  to  the  energy  of  recovery 

d_f  Index  corresponding  to  the  energy  of  muscle  de¬ 

crease  force 

f  Index  corresponding  to  the  energy  of  fatigue 
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Example 

load  lt_splenius 

%  to  compute  the  instantaneous  frequency  and  the 
%  instantaneous  amplitude  slopes 
Fs  =  l 02  4 ; 

sample_interval=60*1024; 

%1  sec  =  1024  samples;  then  60*1024  =  1  minute 
shift=1024;  %shift  every  second 

[p] =linear_regression (fm,  Fs,  sample_interval ,  shift); 
lr_f=p ( : , 1 ) ; 

[p] =linear_regression (i_amp,  Fs,  sample_interval,  shift) 

lr_ia=p (  :  ,  1 )  ; 

[ increase_f orce,  recovery,  decrease_f orce,  fatigue] =... 
jasa_analysis (lr_i,  lr_f) ; 

%  Output  of  the  code  above 

increase_f orce  = 

0.8333 

recovery  = 

1 .6667 

decrease_f orce  = 

62 .5000 
fatigue  = 

35 
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Appendix  B 
Questionnaire 
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Questionnaire  to  report  the  perceived  levels  of  discomfort  reported  by  the  subject  (source  — 
AFRL) 

Q1  Have  you  performed  any  of  the  stretching  techniques  in  the  last  two  hours  ?  (Yes=0,  No=6) 

Q2  According  to  the  scale,  select  a  number  that  corresponds  to  your  physical  state.  (0=No 
Discomfort,  2=Discomfort,  4=Soreness,  6=Pain) 

a  Head 

b  Upper  Neck 

c  Lower  Neck 

d  Shoulders 

e  Upper  Back 

f  Lower  Back 

Q3  According  to  the  scale,  select  a  number  that  corresponds  to  hot  spots  or  accumulation  of 
perspiration  (No  Hot  Spots,  Moderate  Hot  Spots,  Severe  Hot  Spots). 

a  Head 

b  Upper  Neck 

c  Lower  Neck 

d  Shoulders 

e  Upper  Back 

f  Lower  Back 

Q4  According  to  the  scale,  select  a  number  that  corresponds  to  any  numbness  or  loss  of 
sensation.  (0=Normal,  3=Abnormal  Sensation,  6=Numbness) 

a  Head 

b  Upper  Neck 

c  Lower  Neck 

d  Shoulders 

e  Upper  Back 

f  Lower  Back 
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Q5  According  to  the  scale,  select  a  number  that  corresponds  to  a  penetrating  ache  in  your  head. 
(0=No  Headache,  3=Minimal  Headache,  6=Severe  Headache) 

Q6  According  to  the  scale,  select  a  number  that  corresponds  to  your  mental  state.  (0=Relaxed, 
3=Slightly  Tense,  6=Restless) 

Q7  According  to  the  scale,  select  a  number  that  corresponds  to  your  mental  frame  of  mind. 
(0=Alert,  3=Tired,  6=Exhausted) 

Q8  According  to  the  scale,  select  a  number  that  corresponds  to  your  concentration  level. 
(0=Attentive,  3=Distracted,  6=Loss  of  Focus) 

Q9  According  to  the  scale,  select  a  number  that  corresponds  to  the  effort  exerted  to  perform  the 
MVC  at  70%.  (0=No  Effort,  3=Minimal  Effort,  6=Maximal  Effort) 
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Appendix  C 
Matlab  Code 
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%  function  [fltana,  indfltana]=  gfb (n_sb) ; 

%  Purpose: 

%  Generate  the  analysis  filterbank.  For  this  toolbox  the  number  of  subbands 
%  n_sb  can  be  16,  32  or  64.  We  can  represent  the  filter  bank  in  matrix  form 
%  by  h [f,n]=fltana (f*Mf+n+l) for  f =0 , 1 , 2 , . . . F-l  and  n=0 , 1 , 2 , . . . Mf-1  where  F 
%  is  the  number  of  subbands  n_sb  and  Mf  is  the  length  of  the  filter  in 
%  subband  f.  The  length  of  the  filter  f  can  be  computed  from  indfltana 
%  using  Mf=indfltana (f+1) -indfltana (f)  for  f>0,  and  indf ltana ( f +1 )  for  f=0 . 
%  The  output  indfltana  is  required  only  when  the  filters  have  different 
%  length.  In  this  project  all  the  filters  used  had  length  equal  to  128 
%  coefficients. 

%  Synopsis: 

%  [fltana,  indfltana] =  gfb (n_sb) ; 

%  input  vector: 

%  n_sb: number  of  subbands 
%  output  vector: 

%  fltana:  column  vector  with  the  coefficients  of  the  analysis  filter  bank 
%  indfltana:  column  vector  with  the  index  of  the  end  of  each  filter 
%  Example: 

%  fltana,  indf ltana] =gfb (32 )  ; 

%  hn=reshape (fltana, 128, 32) ; 

%  Hw=f ft (hn, 512) ; 

%  w— (-1:2/511:1) ; 

%  plot (w, abs (Hw) ) ;  xlabel ( ' $\f rac { \ omega } { \pi } $ ' , ' interpreter' ,' latex' ) ; 

%  ylabel ('$ | H (\omega) | $',' interpreter' ,' latex' ) ;  grid  on 
%  #Author:  Cristhian  Potes 
%  #Date:  May  23,  2007 

function  [fltana,  indfltana] =  gfb(n_sb); 
if  n_sb==16  |  n_sb==32  |  n_sb==64 

f ilter { 1 , 1 }= [ ' f ilterbank_'  num2str (n_sb)  'subbands']; 
load(filter{l,  1})  ; 
f ilterbank=f ilterbank'  ; 

[N,M] =size (filterbank)  ; 
fltana=f ilterbank (  :  )  ; 
indf ltana= (N : N : M*N) ' ; 
else 

disp ('Number  of  subbands  has  to  be  16,  32,  or  64'); 
f ltana= [ ] ; 
indf ltana= [ ]  ; 
end 
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%  function  [sO, isO] =f_f ilter_bank_ana (sig, fltana, indfltana) 

%  Purpose: 

%  Compute  the  subband  coefficients  sO  by  convolving  the  signal  sig  with  the 
%  filter  coefficients  fltana.  Representing  the  input  signal  as  x[n],  for 
%  n=0, 1, 2, . . . , L-l  where  L  is  the  length  of  the  signal,  then  we  have 
%  x [n] =sig  (n+1)  .  The  output  y[f,n],for  f=0, 1, 2,  .  . . , F-l  and 
%  n=0,l,2, . . . , L+Mf-2  where  F  is  the  number  of  subbands  n_sb  and  Mf  is  the 
%  length  of  the  filter  in  subband  f  is  given  by  y [f , n] =s0 (isO (f ) +n+l)  , 

%  for  f>0  and  sO(n+l),  for  f=0. 

%  Synopsis: 

%  [sO, isO] =f_f ilter_bank_ana (sig, fltana, indfltana) 

%  input  vector: 

%  sig:  column  vector  of  size  LX1  with  the  coefficients  of  the  input  signal 
%  fltana:  column  vector  of  size  F*MfXl  with  the  coefficients  of  the 
%  analysis  filter 
%  bank 

%  indfltana:  column  vector  of  size  FX1  with  the  index  of  the  end  of  each 
%  filter 
%  output  vector: 

%  sO:  column  vector  of  size  sum(L+Mf-l)  for  0<=f<=F+l  with  the  coefficients 
%  of  each  subband 

%  isO:  column  vector  of  size  FX1  with  the  index  of  the  end  of  each  subband 
%  Example: 

%  load  ' NFFO  0  0 1_H1_7  0_EMG' 

%  [sig] =data ( : , 2 ) ; 

%  [fltana,  indfltana] =gfb (32) ; 

%  [sO,  isO] =f_filter_bank_ana (sig, fltana, indfltana) ; 

%  #Author:  Ricardo  von  Borries,  Cristhian  Potes 
%  #Date:  November  13,  2007 

function  [sO, isO] =f_f ilter_bank_ana (sig, fltana, indfltana) 

sig=sig ( : ) ; 

down_samp=l ; 

indsig=length (sig) ; 

indsig= [ 0 ; indsig]  ; 

indfltana= [0 ;  indfltana]  ; 

lisig=length (indsig)  ; 

i0=indsig (lisig-1)  ; 

il=indsig (lisig)  ; 

x=sig (iO+1 : il ) ; 

sig=sig (1 : iO)  ; 

indsig=indsig (1: (lisig-1) ) ; 

lisig=lisig-l; 

for  i=0 : (length (indfltana) -2) 

f=conv (x, fltana ( (indfltana (i  +  1) +1)  : indfltana (i  +  2)  )  )  ; 
f=f ( 1 : down_samp : length ( f )  )  ; 
sig= [sig; f ] ; 

indsig= [indsig; indsig (lisig) flength (f )  ]  ; 
lisig=lisig+l; 
end 

indsig=indsig (2 : length (indsig) )  ; 

isO=indsig; 

s0=sig; 
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%  function  [tfr] =tsub (sO,  isO,  shift); 

%  Purpose: 

%  Arrange  in  a  time  subband  matrix  representation  all  the  subband 
%  coefficients.  This  function  assumes  that  all  the  filters  have  equal 
%  length.  The  time  subband  coefficients  y[f,n]  for  f=0, 1,2, . . . , F-l  and 
%  n=0, 1 , 2, . . . , L+Mf-2  where  F  is  the  number  of  subbands  n_sb  and  Mf  is  the 
%  length  of  the  filter  in  subband  f  is  given  by  y [f , n] =tfr (f+1, n+1 . 

%  Synopsis: 

%  [tfr] =tsub  (sO,  isO,  shift); 

%  input  vector: 

%  sO:  column  vector  of  size  sum(L+Mf-l)  for  0<=f<=F+l  with  the  coefficients 
%  of  each  subband 

%  isO:  column  vector  of  size  FX1  with  the  index  of  the  end  of  each  subband 
%  shift:  vector  of  length  two  with  the  number  of  coefficients  to  remove 
%  from  the  output  sO  to  compensate  the  increase  in  the  length  of  the 
%  subbands  by  the  convolution  operation 
%  output  vector: 

%  tfr:  time  subband  matrix  representation 
%  Example: 

%  load  ' NFFO  0  0 1_H1_7  0_EMG' 

%  [ sig] =data ( : , 2 )  ; 

%  [fltana,  indfltana] =gfb  (32) ; 

%  N=length (fltana) /length (indfltana)  ; 

%  shift= [round (N/2 ), round (N/2 )]  ; 

%  [sO,  isO] =f_filter_bank_ana (sig,  fltana,  indfltana)  ; 

%  [tfr] =tsub (sO,  isO,  shift); 

%  #Author:  Cristhian  Potes 
%  #Date:  July  10,  2007 

function  [tfr] =tsub (sO,  isO,  shift); 
s0=s0 ( : ) ; 
is0=is0  (  :  ) ; 

M=length ( isO ) ; 

L=is0 (1) ; 

tf r=re shape ( sO , L, M) ' ; 
tf r=tf r . ~2 ; 

tfr=tfr (:, shift (1) :L-shift) ;  %  Compensation  of  the  delay 
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%  function  [fm] =inst_freq (tfr) ; 

%  Purpose: 

%  Compute  the  instantaneous  frequency  fm  of  a  given  signal  from  the  time 
%  subband  matrix  representation  tfr.  The  instantaneous  frequency  f[n]  at 
%  time  n=0 , 1 , 2 , . . . , L-l  is  given  by  f [n] =fm (n+1 ) . 

%  Synopsis: 

%  [fm] =inst_freq (tfr)  ; 

%  input  vector: 

%  tfr:  time  subband  matrix  representation 
%  output  vector: 

%  fm:  column  vector  of  size  LX1  with  the  instantaneous  frequency  components 
%  Example: 

%  load  ' NFFO  0  0 1_H1_7  0_EMG' 

%  [ sig] =data (  :  ,  2 )  ; 

%  [fltana,  indfltana] =gfb  (32) ; 

%  N=length (fltana) /length (indfltana) ; 

%  shift= [round (N/2 ), round (N/2 )] ; 

%  [sO,  isO] =f_f ilter_bank_ana (sig,  fltana,  indfltana)  ; 

%  [tfr] =tsub (sO,  isO,  shift); 

%  [ fm] =inst_f req (tsub)  ; 

%  t=  (1 : length (sig) ) /1024;  fm=fm*1024;  plot(t,  fm)  ; 

%  axis([l  max(t)  0  511]);  grid  on; 

%  xlabel ( ' Time [ s ] ' ,  'interpreter',  'latex'); 

%  ylabel (' Frequency [Hz ]' ,  'interpreter',  'latex'); 

%  #Author:  Cristhian  Potes 
%  #Date:  Mayo  22,  2007 

function  [fm] =inst_freq (tfr) ; 

[N, M] =size (tfr)  ; 
freqs= (0 :N-1) *1/ (2*N)  ; 

E=sum (tfr) ; 

fm= (freqs*tfr) ./E; 

fm=fm ( : ) ; 
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%  function  [ i_amp] =inst_amp (tfr ) ; 

%  Purpose: 

%  Compute  the  instantaneous  amplitude  i_amp  of  a  given  signal  from  the 
%  time  subband  matrix  representation  tfr.  The  instantaneous  amplitude  ia[n] 
%  at  time  n=0 , 1 , 2 , . . . , L-l  is  given  by  ia [n] =i_amp (n+1 ) . 

%  Synopsis: 

%  [ i_amp] =inst_amp (tf r ) ; 

%  input  vector: 

%  tfr:  time  subband  matrix  representation  of  size  FXL 
%  output  vector: 

%  i_amp :  column  vector  of  size  LX1  with  the  instantaneous  amplitude 
%  components 
%  Example: 

%  load  ' NFFO  0  0 1_H1_7  0_EMG' 

%  [ sig] =data (  :  ,  2 )  ; 

%  [fltana,  indfltana] =gfb (32) ; 

%  N=length (fltana) /length (indfltana)  ; 

%  shift= [round (N/2 ), round (N/2 )]  ; 

%  [sO,  isO] =f_filter_bank_ana (sig,  fltana,  indfltana)  ; 

%  [tfr] =tsub (sO,  isO,  shift); 

%  [ i_amp] =inst_amp (tsub) ; 

%  t= (1 : length (sig) ) /1024;  plot (t,  i_amp) ; 

%  axis([l  max(t)  0  max ( i_amp) ] ) ;  grid  on; 

%  xlabel ( ' Time [ s ] ' ,  'interpreter',  'latex'); 

%  ylabel (' Amplitude' ,  'interpreter',  'latex'); 

%  #Author:  Cristhian  Potes 
%  #Date:  Mayo  25,  2007 

function  [ i_amp] =inst_amp (tf r) ; 
ws=l ; 
s  s  =  1 ; 

[N, M] =size (tfr)  ; 
ia=zeros (M, 1 ) ; 
i_amp=ia; 
for  1=1 : N 

nc=ceil (2*N/ss) ; 
ss=ss+l; 

[is,  ie] =windowing (nc,  M,  ws ,  0 )  ; 
ana_sig=tfr (i,  :  )  ; 
for  j=l : length (is) 
ww=abs (ana_sig (is ( j ) : ie ( j ) ) ) ; 
ia  (is ( j ) ) =max (ww)  ; 
end 

ia=ia ( : ) ; 
i_amp= i_amp + i a ; 
ia=zeros (M, 1 ) ; 

end 

i_amp=sqrt (i_amp) ; 
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%  function  [ index_start ,  index_end] =windowing (nc,  L,  , shift,  option) 

%  Purpose: 

%  Compute  the  start  and  end  indexes  of  the  location  of  a  rectangular  window 
%  for  a  specified  shift  amount.  This  function  is  useful  when  we  want  to  apply 
%  a  rectangular  window  to  a  signal  of  length  L 
%  Synopsis: 

%  [ index_start ,  index_end] =windowing (nc,  L,  shift,  option);  %option  default  =  1 
%  [ index_start ,  index_end] =windowing (nc,  L,  shift);  %shift  default  =1 
%  [ index_start ,  index_end] =windowing (nc,  L) ; 

%  input  vector: 

%  nc:  size  of  the  rectangular  window,  for  l<=nc<=L 
%  L:  length  of  the  signal 

%  option:  if  1  ->  ( index_end-index_start )  is  always  nc.  The  rectangular 
%  window  is  not  zero  padded  at  the  end  of  the  signal 

%  if  0  ->  ( index_end-index_start )  is  not  always  nc.  The  rectangular 

%  window  is  zero  padded  at  the  end  of  the  signal 
%  shift:  shifted  window,  for  l<=shift<=nc; 

%  output  vector: 

%  index_start :  column  vector  with  the  start  indexes; 

%  index_end:  column  vector  with  the  end  indexes; 

%  Example: 

%  %  Let  us  suppose  that  we  want  to  apply  a  rectangular  window  of  5  points  with 
%  a  shift  of  3  to  a  signal  of  length  15 
%  nc=5;  L=15;  shift=3; 

%  when  option=l 
%  option=l; 

%  [ index_start , index_end] =windowing (nc,  L,  option,  shift) 

%  when  option=0 
%  option=0; 

%  [ index_start , index_end] =windowing (nc,  L,  option,  shift) 

%  #Author:  Cristhian  Potes 
%  #Date:  Feb  16,  2007 

function  [index_start,  index_end] =windowing (nc,  L,  shift,  option) 

if  nargin==3 
option=l ; 
end 

if  nargin==2 
shift=nc; 
option=l ; 
end 

if  nc>L 

disp(' Length  of  the  rectangular  window  nc  has  to  be  l<=nc<=L' ) ; 
end 

if  shift>nc  |  shift==0 

disp(' Shift  of  the  window  has  to  be  l<=shift<=nc'  )  ; 
end 

index_start=l : shift : L; 
index_end=nc : shift : L; 

dif=abs (length (index_end) -length ( index_start ) )  ; 
if  dif  ~=0 
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index_end= [ index_end,  ones (1,  dif).*L]; 

end 

if  option  ==1 

c=find (index_end-index_start+l==nc) ; 
index_start=index_start (c) ; 
index_end=index_end (c)  ; 
end 

index_start=index_start ( : ) ; 
index_end=index_end ( : )  ; 
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%  function  [p] =linear_regression (s,  Fs,  ns,  shift) 

%  Purpose: 

%  Estimate  the  slopes  (a_i)  and  points  of  interception  (b_i)  of  a  piecewise 
%  approximation  a_i*t+b_i,  for  the  time  interval  t_iO<=t<=t_il  of  the  input 
%  data. 

%  Synopsis: 

%  [p] =linear_regression (s,  Fs,  ns,  shift);  shift  default=l; 

%  [p] =linear_regression (s,  Fs,  ns); 

%  input  vector: 

%  s : input  data 
%  Fs:  sampling  frequency 

%  ns:  size  of  the  interval  t_iO<=t<=t_il 

%  shift:  shifted  window,  for  l<=shift<=length ( s ) ;  The  intervals  of  the  linear 
%  approximation  may  overlap 
%  output  vector: 

%  p:  matrix  where  the  first  column  corresponds  to  the  slopes  and  the  second 
%  column  to  the  points  of  interception 
%  Example: 

%  load  lt_splenius  to  compute  the  instantaneous  frequency  slopes 

%  s=fm;  %  the  instantaneous  frequency  fm  is  the  output  of  the  function  inst_freq.m 
%  Fs=1024 ; 

%  ns=60*1024;  %1  sec  =  1024  samples;  then  60*1024  =  1  minute 
%  shift=1024;  %shift  every  minute 
%  [p] =linear_regression (s,  Fs,  ns,  shift); 

%  t=round ( int_time/2 ) + (0 : size (p, 1) -1) *shift; 

%  t=t . / (60*1024) ; 

%  plot(t,  p ( : , 1 ) , ' r . ' , ' MarkerSize' , 15 ) ;  grid  on 

%  xlabel (' Time [min] interpreter' ,  'latex',  ' f ontweight ' ,  'bold',... 

%  ' font size' ,  11) ; 

%  ylabel (' Slope (Hz/min) ' ,  'interpreter',  'latex',  ' f ontweight ' ,  'bold',... 

%  ' font size' ,  11) ; 

%  #Author:  Cristhian  Potes 
%  #Date:  Aug  31,  2007 

function  [p] =linear_regression (s,  Fs,  ns,  shift) 

if  nargin==3 
shift=ns ; 
end 

s  =  s  (  :  )  ; 

N  =  length ( s ) ; 

time= ( 1 : N) . /Fs ;  time=time ( : ) ; 

[n,  m] =windowing (ns ,  N,  shift); 

for  i  =  1: length (n) 
x  =  s (n  (i)  :m (i)  )  ; 
t  =  time (n (i) :m (i) ) ; 

P  ( i ,  : )  =  polyfit (t,  x,  1); 

end 
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%  function  [ increase_force,  recovery,  decrease_force,  fatigue] = 

%  jasa_analysis (s_inst_amp,  s_inst_freq) 

%  Purpose: 

%  computes  a  joint  analysis  of  the  spectrum  and  the  amplitude  of  the 
%  electromyography  signals  ( JASA) from  the  instantaneous  frequency  and  the 
%  instantaneous  amplitude  slopes  to  compute  indexes  of  muscle  increase  force, 

%  recovery,  muscle  decrease  force  and  fatigue. 

%  Synopsis: 

%  [ increase_f orce,  recovery,  decrease_force,  fatigue] = 

%  jasa_analysis (s_inst_amp,  s_inst_freq) 

%  input  vector: 

%  s_inst_amp:  instantaneous  amplitude  slopes 
%  s_inst_freq:  instantaneous  frequency  slopes 
%  output  vector: 

%  if:  index  corresponding  to  the  energy  of  muscle  increase  force 
%  r:  index  corresponding  to  the  energy  of  recovery 
%  df :  index  corresponding  to  the  energy  of  muscle  decrease  force 
%  f:  index  corresponding  to  the  energy  of  fatigue 
%  Example: 

%  load  lt_splenius  to  compute  the  instantaneous  frequency  and  the 
%  instantaneous  amplitude  slopes 
%  Fs=1024 ; 

%  sample_interval=60*1024;  %1  sec  =  1024  samples;  then  60*1024  =  1  minute 
%  shift=1024;  %shift  every  minute 

%  [p] =linear_regression (fm,  Fs,  sample_interval ,  shift);  lr_fm=p ( : , 1 ) ; 

%  [p] =linear_regression (i_amp,  Fs,  sample_interval ,  shift);  lr_inst_am [=p ( : , 1 ) ; 

%  [ increase_f orce,  recovery,  decrease_force,  fatigue]  = 

%  jasa_analysis (lr_inst_amp,  lr_fm) 

%  #Author:  Cristhian  Potes 
%  #Date:  Nov  13,  2007 

%  #Date  modified:  March  10,  2008 

function  [increase_force,  recovery,  decrease_force,  fatigue]=... 
jasa_analysis (s_inst_amp,  s_inst_freq) 

if  length ( s_inst_freq) >0  &  length ( s_inst_amp) >0 

if_ia=f ind ( s_inst_amp>0 ) ; 

increase_f orce_ia=s_inst_amp (if_ia) ; 

increase_f orce_if=s_inst_f req (if_ia) ; 

if_if =f ind ( increase_f orce_if >0 ) ; 

increase_f orce_if=increase_f orce_if (if_if) ; 

increase_f orce_ia=increase_f orce_ia (if_if ) ; 

increase_f orce=increase_f orce_ia . " 2+increase_f orce_if . "2 ; 

increase_f orce=sum (increase_force) ; 

if_ia=f ind ( s_inst_amp<0 ) ; 
recovery_ia=s_inst_amp ( if_ia) ; 
recovery_if=s_inst_f req ( if_ia) ; 
if_if=f ind (recovery_if>0 ) ; 
recovery_if =recovery_if (if_if) ; 
recovery_ia=recovery_ia (if_if ) ; 
recovery=recovery_ia  .  ,'2+recovery_if  .  "2; 
recovery=sum (recovery)  ; 
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if_ia=f ind ( s_inst_amp<0 ) ; 

decrease_f orce_ia=s_inst_amp (if_ia) ; 

decrease_f orce_if=s_inst_f req (if_ia) ; 

if_if=f ind (decrease_f orce_if <0 ) ; 

decrease_f orce_if=decrease_f orce_if (if_if) ; 

decrease_f orce_ia=decrease_f orce_ia  (,i'f_if )  ; 

decrease_f orce=decrease_f orce_ia . " 2+decrease_f orce_if . ~2 ; 

decrease_f orce=sum (decrease_force)  ; 

if_ia=f ind ( s_inst_amp>0 ) ; 
f atigue_ia=s_inst_amp (if_ia) ; 
f atigue_if=s_inst_f req ( if_ia)  ; 
if_if =f ind ( f atigue_if <0 ) ; 
f atigue_if=f atigue_if (if_if ) ; 
f atigue_ia=f atigue_ia (if_if ) ; 
f atigue=f atigue_ia . "2  +  f atigue_if . " 2 ; 
fatigue=sum (fatigue)  ; 

else 

increase_f orce=0 ; 
recovery=0 ; 
decrease_f orce=0 ; 
f atigue=0 ; 
end 
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